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SECTION  1 


INTRODUCTION 

This  document  is  the  final  report  for  Contract 
Number  N00014 - 81 -C - 2041  ,  which  covered  work  performed  for 
Code  4707  of  the  Naval  Research  Laboratory  (NRL)  by  the 
Plasma  Physics  Division  of  Science  Applications,  Inc.  (SAI) 
during  the  period  24  November  1980  to  24  September  1982. 

The  material  covered  in  this  report  consists  of 
three  general  areas  in  which  plasma  physics  plays  a  signifi 
cant  role  in  the  modeling  of  radiation  sources  for  the 
advanced  simulation  research  program  sponsored  by  the 
Defense  Nuclear  Agency  (DNA) .  The  first  is  the  description 
of  a  basic  model  for  the  implosion  of  a  system  of  identical 
wires  driven  by  a  pulsed-power  generator.  The  second  is  a 
model  for  computing  the  linear  ideal  MHD  instability  growth 
rates  for  azimuthally- symmetric ,  cylindrical  Z-pinch  equili 
bria.  This  analysis  includes  both  kink  and  sausage-type 
perturbations  of  the  equilibrium.  The  third  area  concerns 
the  properties  of  magnetically-insulated  power  feeds  for 
driving  imploding  Z-pinch  loads. 


SECTION  2 


THE  "WIRES"  CODE  --  A  SIMPLE  MODEL  FOR  IMPLOSIONS 

It  has  been  known  for  several  years  that  imploding 
wire  arrays  during  run-in  obey  the  so-called  "F-ma"  dynamics 
to  remarkable  accuracy.1  Recent  data  has  suggested  that  for 
truly  massive  arrays  (~  500  ug)  the  F  =  ma  scaling  finally 
breaks  down,  with  arrays  apparently  going  unstable  prior 
to  achieving  a  significant  inward  acceleration.  For  arrays 
of  <  300  ug,  however,  the  F  =  ma  formalism  appears  good, 

2 

and  predicts  the  assembly  time  to  within  a  few  nanoseconds. 
Given  that  the  dynamics  of  the  wire  centers  has  been  under¬ 
stood  at  this  level,  it  appears  reasonable  to  attempt  a 
generalization  of  the  F  =  ma  model  to  include  a  radiation 
algorithm  and  a  prescription  for  tracking  the  internal 
energy  of  the  wires. 

One  strong  motivation  for  pursuing  this  type  of 
model  for  the  early-time  behavior  of  wire  arrays  is  to 
establish  the  correct  initial  conditions  for  one-dimen¬ 
sional,  radiation-coupled  hydro  codes,  such  as  WHYRAD  and 
SPLATT,  which  currently  assume  that  the  individual  wires 
instantly  expand  into  a  plasma  annulus  with  a  temperature  of 
1-10  eV  prior  to  implosion.  By  calculating  the  intial 
conditions  from  the  generalized  F  =  ma  model  described 
here,  these  codes  will  be  able  to  provide  scaling  with 
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number  of  wires  varied  and  total  array  mass  fixed,  which 
currently  is  not  possible  to  compute.  Also,  the  results 
of  many  runs  of  h'HYRAD  and  SPLATT  indicate  that  during 
this  early  implosion  period,  the  radiation  is  basically  a 
black-body  spectrum,  and  the  plasma  remains  relatively 
cool,  indicating  that  the  detailed  radiation  and  chemistry 
package  of  these  sophisticated  codes  is  not  needed  for  the 
run-in  phase  of  the  implosion.  By  using  the  simple  model  de 
scribed  here,  the  early-time  behavior  of  the  array  can  be 
obtained  in  a  small  fraction  of  the  computer  time  required 
in  KHYRAD,  thereby  allowing  the  detailed  models  to  be 
utilized  where  they  are  most  necessary,  during  the  assembly 
and  compression  of  the  plasma  annulus  on  axis. 

Figure  2-1  shows  a  schematic  representation  of  the 
wire  array  at  timet.  The  individual  wires  have  radius, 
a(t),  the  wire  array  has  radius,  r(t),  and  the  entire 
system  of  N  wires  is  enclosed  by  a  cylinder  of  radius,  b, 
which  carries  the  return  current.  The  external  circuit  is 
shown  in  Figure  2-2,  and  consists  of  an  external  voltage 
generator,  providing  a  voltage  waveform,  V(t),  with  a 
generator  impedance  ZQ.  This  generator  section  could  be 
replaced  by  a  transmission  line  of  impedance,  ZQ  and  length, 
t,  which  is  initially  fully  charged.  The  generator  drives 
a  time-dependent  load  described  by  the  diode-housing  induc¬ 
tance,  Lp,  and  the  time-varying  plasma  resistance  and 
inductance,  R  ( t )  and  L  (t)  respectively. 
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The  plasma  circuit  parameters  are  assumed  to  be 
correctly  given  by  the  Russel  formula  for  inductance  and 
the  Spitzer  resistivity, 


for  lengths  in  centimeters ,  where  l  is  the  array  length, 
and 
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is  the  Spitzer  resistivity  for  electron  temperature,  T  in 
Kelvins,  ^  A  being  the  Coulomb  logarithm  and  y£  a  factor  of 
order  unity  which  depends  only  on  th.e  effective  ionization 
state,  Zef£.  For  an  element  of  atomic  number,  Z,  the 
effective  ionization  state  is  approximately, 
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and  the  current  flowing  in  each  wire  is  I  (t)/N. 

The  motion  of  the  array  is  given  by  the  J  x  B  force 
for  each  wire.  If  each  wire  has  mass/length  *  y,  the 
radius  of  the  array  is  given  by 
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which  is  integrated  numerically  as  two  first-order  equations 
together  with  the  circuit  equation,  using  a  Runge-Kutta 
integrator . 

The  radiated  power  is  modeled  as  a  black-body  with 
emissivity  given  by 

>,->  <r< 

e  =  e  f  +  e  f  , 
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where  e  (c  )  is  the  emissivity  for  photons  greater  (less) 
than  a  specified  cut-off  energy,  E*,  and  f^"  (f<)  is  the 
fraction  of  the  blackbody  output  which  is  emitted  above 
(below)  E*.  Clearly,  f >  +  f <  =  1 .  The  total  radiated  energ} 
"rad’  is  t^ien  given  by 


dw  , 
rad 

~~xr~ 

and  the  yield, 

dw>  , 
rad 


=  eoT^A 

s 

"rad’  ab°ve  E*  is  given  by 
r>  >  ~4. 

*  f  c  ffT  A  , 


where  a  is  the  Stefan-Boltzmann  constant  and  A„  =  2 tt a fi. N  is 

s 

the  total  surface  area  of  the  array. 

To  complete  the  description  of  the  model,  a  prescrip¬ 
tion  for  determining  the  plasma  temperature,  T,  and  the 
wire  radius,  a,  is  required.  If  each  wire  is  assumed  to  be 
a  Bennett  equilibruim,  the  Bennett  pinch  condition, 


nd  -  Zeff)  KbT 


9 


provides  a  relationship  between  temperature  and  current, 


is  Boltrmanr. ’ s  constant. 
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where  M  is  the  atomic  mass  and  Kg 
The  wire  radius,  a,  may  be  obtained  from  the  energy  balance 
between  Ohmic  dissipation  and  radiation, 

I^R  =  eaT4A 
P  p  .  s 


or 
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The  model  described  above  can  be  integrated  until 
the  wires  iust  touch,  a  =  r  sin(r/.\),  at  which  point  the 
system  of  individual  wires  coalesces  into  a  plasma  annulus, 
which  rapidly  assembles  on  axis  converting  the  kinetic 
energy  of  implosion  to  temperature,  radiation  and  out¬ 
going  kinetic  energy. 

For  simple  scaling  law  studies,  the  following  very 
crude  model  has  been  implemented  to  model  the  assembly. 

The  plasma  annulus  is  converted  to  a  cylinder  of  radius, 
r  ,  as  shown  in  Figure  2-5,  with 
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where  a  is  the  wire  radius  when  the  wires  just  touch. 

The  plasma  temperature  is  adjusted  so  that  the  kinetic  energy 
is  entirely  absorbed  into  temperature, 

(1  *  zeff)  AT  -  i  N,vJ/Kb  . 

The  system  is  then  allowed  to  radiate  and  cool  for  a  period 
of  five  MHD  growth  times,  calculated  as  Alfven  transit  times, 

t  -  S(ro/vA). 

2  k 

where  V.  =  (B  /4-jp)^  is  the  Alfven  speed,  B  is  the  magnetic 

field  at  r  =  rQ  due  to  the  current  I  and  p  is  the  mass 

2 

density,  p  =  Ny/iir  .  During  the  cooling  period,  assuming 
blackbody  radiation,  the  code  separately  integrates  for 
the  total  yield  and  the  yield  above  E*. 

A  test  case  has  been  run  for  an  Aluminum  (1=15, 
A=27)  array  driven  at  constant  voltage,  V(t)=V  ,  with 
the  following  parameters: 

N  *  Number  of  Wires  =  6 
Ny£  =  Array  Mass  =  100  yg 
£  =  Array  Length  =  3  cm 
r(o)  =  Initial  Array  Radius  =  2.2  cm 
b  =  Return-Current  Radius  =  5  cm 
V  =  Open-Circuit  Voltage  =  3  MV 
ZQ  *  Generator  Impedance  =  0.7  n 


=  Diode -Hous ing  Inductance  =  10  nH 

e>  =  Emissivity  for  hv>lkev  =  5  x  10  6 

<  -  4 

c  =  Emissivity  for  hv<lkev  =  5  x  10 

The  characteristics  of  the  implosion  are  summarized  in 
Figures  2-4  thru  2 -7.  As  seen  in  Figure  2-4,  the  implo¬ 
sion  of  this  array  requires  approximately  69  ns.  At  the 
time  the  wires  touch,  as  in  Figure  2-3,  the  wires  have 

g 

achieved  an  inward  speed  of  1.3  x  10  cm/sec.  The  individ¬ 
ual  wire  radius,  a,  varies  over  a  factor  of  two  during 
most  of  the  implosion.  At  very  early  times,  this  radius 
is  artifically  large  due  to  the  assumptions  of  constant 
emissivities  and  the  prescription  of  choosing  "a"  as  the 

radius  where  Ohmic  heating  is  balanced  by  radiative  cooling. 

2 

Experiments  at  Maxwell  have  displayed  an  initial  pinching 
of  the  individual  wires  followed  by  an  expansion  of  the 
wires,  which  is  at  least  qualitatively  as  shown  in  these 
calculations . 

Figure  2-5  shows  the  temperature  and  average 
ionization  state  vs.  time  for  this  implosion.  The  tempera¬ 
ture,  which  is  tied  to  the  current  by  the  Bennett  pinch 
condition,  peaks  at  approximately  655  eV  at  peak  current, 
and  subsequently  drops  to  249  eV  by  the  end  of  the  run-in. 
The  ionization  state,  Zeff>  is  between  9  and  11  during 
most  of  the  implosion. 


a 
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equation  by  multiplying  both  sides  by  the  total  current. 

I  ,  to  obtain 
P 


V  =  3T  [t  <ld*lp!Ip]  * 


(Z  +R  'll" 

o  p-  p 


~  l2l 

2  p  p 


where  IV  is  the  input  power  from  the  external  generator, 

which  must  equal  the  rate  at  which  energy  is  stored  in  the 

magnetic  field,  Ohmic  power  losses,  and  the  rate  at  which 

energy  is  stored  as  kinetic  energy  of  the  array  (the 

term).  Figure  2-6  illustrates  how  these  various  components 

of  the  power  equation  vary  in  time.  The  rate  at  which 

energy  is  stored  in  the  magnetic  field  is  not  plotted,  but 

is  iust  the  difference  between  the  V  I  curve  and  the  sum 

o  p 

of  the  other  two  curves.  At  the  end  of  the  run-in,  the 
wires  are  acquiring  kinetic  energy  at  a  rate  which  exceeds 
VQI  ,  and  in  fact  the  implosion  is  tapping  stored  field 
energy  iust  prior  to  assembly  on  axis. 

Figure  2-"  shows  the  evolution  of  the  various  energy 
channels  during  the  implosion.  During  run-in,  when  the 
temperature  is  low,  internal  energy  and  radiation  are  rela¬ 
tively  small  compared  with  field  energy  and  kinetic  energy. 
Again,  at  the  end  of  the  run-in,  the  field  energy  decreases 
rapidly  as  the  kinetic  energy  increases. 


After  the  wires  touch,  at  t  =  69  ns,  the  code 

instantly  converts  the  annulus  to  a  cylinder  and  "shock 

heats"  it  by  converting  all  the  kinetic  energy  to  internal 

energy.  This  prescription  in  the  test  case  yields  a  cylinder 

°0  -  i 

of  0."  cm  diameter  with  an  ion  density  of  2  x  10“  cm  "  at  a 
temperature  of  11.4  keV.  The  cylinder  cools  rapidly  by 
radiative  power. loss,  and  after  five  Alfven  transit  times 
(or  5.1  ns)  its  temperature  has  dropped  to  113  eV.  In  this 
problem,  the  total  energy  radiated  at  all  frequencies  is 
34.2  kJ  and  the  energy  above  1  keV  is  6.7  kJ ,  assuming  a 
blackbody  spectrum.  During  the  run-in,  the  individual 
wires  radiated  28.2  kJ  at  all  frequencies,  but  only  1.7  kJ 
above  1  keV.  The  radiation  above  1  keV,  therefore,  occurs 
after  collapse  in  this  model,  but  a  significant  fraction 
of  the  total  yield  can  occur  during  the  run-in.  The  model 
assumes,  of  course,  that  the  individual  wires  remain  stable 
and  do  not  develop  "hot  spots"  during  the  collapse.  If  hot 
spots  develop,  they  may  cause  a  larger  fraction  of  the  yield 
above  1  keV  to  occur  during  the  run-in  than  is  calculated 
here . 

Treating  the  test  case  described  above  as  a  base 
case,  a  parameter  study  has  been  made  to  test  the  sensi¬ 
tivity  of  the  radiation  yield  to  variations  of  several  of  the 
input  parameters.  The  results  of  this  study  are  shown  in 
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Table  2-1  and  in  Figures  2-3  thru  2-15.  Multiple  para¬ 
meter  variations  from  the  base  case  have  not  been  attempted, 
but  rather  only  a  single  parameter  has  been  altered  for  each 
of  these  runs.  In  each  case  the  total  radiation  yield, 
and  the  yield  h’* ,  for  hv  > 1  keV,  is  presented.  These  are 
calculated  as  described  above,  using  a  blackbody  spectrum 
with  different  emissivities  above  and  beiow  1  keV,  and 
including  the  radiation  from  a  collapsed,  "shock-heated" 
plasma  cylinder  as  it  cools  during  five  Alfven  transit  times 

Figure  2-8  shows  the  variation  with  X,  keeping  the 
total  array  mass  fixed  at  100  yg.  As  the  number  of  wires 
increases,  the  time  for  impact  is  shortened,  thereby  reduc¬ 
ing  the  total  yield  by  more  than  two-fold.  The  yield  above 
1  keV,  however,  becomes  very  insensitive  to-  the  number 
of  wires  for  N>12. 

Figure  2-9  shows  the  effect  of  varying  the  total 
array  mass,  with  the  number  of  wires  fixed  at  N=6.  Here, 
the  heavier  arrays  are  worse  as  expected  in  this  model. 
Increasing  the  initial  array  radius,  r(o),  as  shown  in 
Figure  2-10,  improves  the  yield.  The  larger  arrays  (at 
100  ug)  take  longer  to  collapse,  acquiring  more  kinetic 
energy.  Also,  since  the  current  return  has  been  fixed  at 
b  =  5  cm,  the  initial  inductance,  L  ,  for  the  largvr  arrays 

P  i 

is  smaller,  thereby  slightly  reducing  the  current  risetime. 


Figures  2-11  thru  2-15  show  the  effect  of  variations 


of  the  external  circuit  parameters .  Increasing  the  open- 

circuit  voltage,  Figure  2-11,  or  reducing  the  generator 

impedance,  Figure  2-12,  both  strongly  improve  the  yield. 

While  this  trend  suggests  that  the  yield  may  simply  depend 

2 

monotonicallv  on  the  power,  VQ/Z  ,  independent  of  whether 

V  or  ZQ  is  the  quantity  varied,  a  detailed  look  at  Table 

2-1  shows  that  this  proposition  is  not  correct.  The  yield 

for  ZQ  =  1.5  Vi  and  VQ  =  3 MY  (V^/Z  =  6  TW)  is  significantly 

lower  than  the  vield  for  2  =  0.7  Cl  and  V  =  2MV  (V'VZ  = 

•  o  o  o  o 

5 . 7  TW) . 

Figure  2-13  shows  that  the  yield  is  insensitive  to 
the  diode -housing  inductance  in  the  range  10  nH  to  20  nH. 

This  insensitivity  is  probably  due  to  the  plasma  inductance, 
which  varies  from  nominally  5  nH  to  50  nH  during  the  implosion 
This  model  has  been  utilised  to  provide  initial  con¬ 
ditions  for  the  SQUEEZE  code,  which  computes  the  collapse  of 
a  plasma  annulus.  The  WIRES  model,  described  above,  is  run 
until  the  individual  wires  iust  touch.  A  bridge  subroutine 
is  then  employed  to  convert  the  wire  array  to  an  imploding 
plasma  annulus.  The  wire  array,  consisting  of  N  wires  of 
radius,  a^,  on  a  circle  of  radius,  v^ ,  is  converted  to  an 
annulus  with  outer  radius,  r  ,  and  inner  radius,  r.,  given  bv. 

'  O  1 
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lhe  SQUEEZE  code  then  continues  the  calculation,  using  a 
more  sophisticated  radiation  and  hydrodynamics  model. 
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Figure  2-3.  Conversion  of  Wire  Array  to  Plasma  Cylinder 
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Figure  2-4.  Array  Radius,  fndi vidua  I -Wire  Radius  and  Implosion 
Snood  vs.  Tinio  for  tho  Rnso-rnsr* 


Level 


Figure  2-6.  Input  Power,  Ohmic  Power,  and  Power  to  Kinetic  Hnergy 
vs.  Time  for  the  Rase  Case 


Radiation 


SECTION  3 


LINEAR,  IDEAL  MHD  STABILITY  ANALYSIS 

Experiments  on  imploding  wire  arrays,  gas  puffs  and 
foils  have  displayed  hot  spots,  beads,  plasma  jets  and  kinks, 
all  of  which  are  believed  to  be  manifestations  of  MHD  in¬ 
stabilities.  These  phenomena  couple  strongly  to  the  plasma 
dynamics  and  may  actually  determine  the  strength  of  the  pinch 
and  the  time  duration  of  the  assembled  plasma.  The  densities 
and  high  temperatures  of  the  hot  spots  or  beads  may  provide 
the  conditions  needed  for  generating  most  of  the  radiation  above 
1  keV. 

The  usual  simple  test  for  the  importance  of  MHD 
instabilities  in  a  plasma  system  is  whether  the  time  required 
for  an  Alfven  wave  to  cross  the  system  is  short  compared  with 
the  confinement  time  of  the  system.  The  Alfven  speed  is  given 
by 

2  I  / 

vA  =  (BV4ttp)A/“ 

where  B  is  the  magnetic  field  and  p  is  the  mass  density.  To 
make  this  argument  specific,  assume  that  a  wire  array  with  a 
mass  of  100  ug  has  collapsed  to  a  plasma  cylinder  of  0.1  cm 
radius  and  3  cm  length,  and  is  carrying  a  current  of  2  MA. 

For  these  parameters,  which  are  typical  of  experimental  con¬ 
ditions,  the  mass  density  is  c  £  10  ‘^g/cm^  and  the  magnetic 
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field  at  the  edge  of  the  cylinder  is  B  *  4  MG,  which  yields 
v^  =  3.6  x  10  cm/sec  and  the  Alfven  transit  time  across  the 
plasma  radius  is  2.8  ns,  which  is  only  about  10°o  of  the 
observed  radiation  pulsewidth  and  the  observed  plasma 
confinement  time.  MHD  instabilities  are  therefore  expected 
to  be  important  in  this  system. 

During  the  implos i on, currents  penetrate  into  the 
plasma  wires  and/or  annulus,  and  the  MHD  growth  rate  will 
be  sensitive  to  the  actual  current  distribution  in  the 
plasma.  The  radiation-coupled  hydro  codes,  WHYRAD  and 
SPLATT,  model  the  field  penetration  and  can  provide  an 
"equilibrium"  configuration  for  the  assembled  plasma.  During 
the  run-in  phase  of-  the  implosion,  inertial  terms  in  the 
zero-order  equations  will  be  important.  A  formulation  for 
the  MHD  instability  growth  rate  for  cylindrical  MHD  equilibria 
is  described  in  this  section  for  an  arbitrary  equilibrium 
current  distribution  which  is  consistent  with  equlibrium 
pressure  balance.  Several  references'^  ^  discuss  this  type 
of  model. 

The  linearized  equations  for  ideal  MHD  may  be  written 
in  terms  of  density  p,  velocity  v,  pressure  p,  current 
density  J,  and  magnetic  field  B  as 


e  Jt  =  Vp  +  ±  x  £  +  ±  x  E 

V  x  B  =  ^  J 
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where  superscripts  "o"  and  "1"  denote  zero-order  and  first- 
order  quantities,  respectively,  and  T  is  the  ratio  of  specifi 
heats . 

These  equations  may  be  expressed  as  a  second-order 
equation  for  the  displacement  vector,  £(x,t),  defined  as 
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This  formula  is  the  starting  place  for  all  linear,  ideal 
MHD  stability  analyses. 

For  a  circular  cylinder  equilibrium,  the  coefficient 
in  the  lineariced  MHD  equations  are  independent  of  6  and  2. 
Each  Fourier  harmonic  of  the  perturbation  will  therefore 
evolve  independently,  and  the  perturbation  may  be  expressed 
as 


C(£,t)  =  i(r)ei(k!"m6'“t) 

In  this  case  it  has  been  shown^  that  the  problem  reduces  to 
a  single,  homogeneous,  second-order  equation  of  the  eigen¬ 
functions  associated  with  the  radial  displacement,  ^ .  This 
equation  is  given  by 


(a£r)  +  qCr  =  0 


where  prime  indicates  differentiation  with  respect  to  the 

radial  coordinate.  The  coefficients,  a  and  q,  are  given  by 

rAC 
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detA  /AllV 
r  |_aca12  \x12J  J 


where  A  is  a  2  x  2  matrix  of  the  form, 


J  ‘‘A  _ 
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For  equilibria  with  no  flow  and  without  an  axial  magnetic 
field,  i.e.  B  =  B_  ,  these  quantities  may  be  expressed  in 
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terms  of  the  Alfven  speed,  \*A,  and  the  sound  speed,  c s , 
where 
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and  p*  denotes  the  total  equilibrium  pressure,  which  satisfies 


The  functions  A.C.A^.A^-,  and  A?1  are  therefore  algebraic 
functions  of  the  eigenvalue  parameter,  w",  and  the  equilibrium 
fields  . 

The  matrix,  A,  contains  the  second  derivative  of  p*, 
and  the  coefficient,  q,  requires  the  derivative  of  A.^/A^. 
These  non-algebraic  features  can  be  troublesome,  particularly 
when  the  equilibrium  data  are  obtained  numerically.  The 
numerical  computation  of  these  derivatives  is  expected  to 
require  some  smoothing  of  the  equilibrium  data. 

The  boundary  conditions  on  £r  are  that  it  vanish  at 
the  outer  boundary,  assumed  to  be  a  conducting  wall,  while 
regularity  at  the  origin  may  be  used  to  determine  its  behavior 
at  r  =  0  from  an  indicial  equation.  Writing  £r  *  rp 
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the  solution  near  the  origin  satisfies  u^=l  for  m=0  and 

?  7 

(u+l)-  =  nr  for  m^O. 

With  this  formulation,  the  problem  is  completely  posed 

in  terms  of  functions  of  the  equilibrium  fields.  In  ideal 

? 

MHD,  the  eigenvalue  is  oT,  implying  solutions  that  are  either 
purely  oscillatory  or  purely  growing,  a  feature  which  follows 
from  the  self-adjoint  nature  of  the  perturbed  fluid  equations. 
More  generally,  eg  ■  for  equilibria  with  flow,  the  self-adjoint 
property  is  lost,  and  the  eigenvalues  will  be  complex. 


The  cylindrical  MUD  stability  problem,  as  formulated 
above,  may  be  solved  numerically  by  a  "shooting  code".-3  In 
this  technique  one  selects  the  equilibrium  and  sets  r  near 

the  origin  to  satisfy  the  indicial  equation  for  a  trial  value 

.  ->  -> 
of  to-.  By  solving  repeatedly  for  £  (r)  for  different  to”, 

a  value  is  found  for  which  5  (wall)  =  0,  the  outer  boundary 
condition.  This  technique  is  a  widely  used  approach. 

In  the  limit  of  surface  currents,  the  linear  stabilit) 
problem  can  be  solved  exactly  and  analytically.  The  eigen¬ 
value  is  given  by 

2 

w2  *  {-£)  [-Sjka  +  m282] 

where  v^  is  the  Alfven  speed  at  the  edge  of  the  plasma  (r=a) 
and  a  perfectly  conducting  wall  is  assumed  at  r=b .  The  co¬ 
efficients,  and  B?,  are  given  by 
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where  I  and  K  are  modified  Bessel  functions, 
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The  ideal  MHD  linear  growth  rate  is  plotted  against 
ka  lor  various  azimuthal  mode  numbers,  m,  and  various  b  a 
values  in  Figure  5-1.  The  results  appear  quite  insensitive 
to  b/a  for  b/a  >  5.  As  b/a  nears  unity,  however,  the  in¬ 
stabilities  with  m>0  are  stabilized  by  wall  stabilization. 

On  the  figure,  an  instability  is  indicated  by  a  negative 
eigenvalue,  i .  e .  <  0.  The  m=0  sausage  mode  is  always 

unstable.  The  modes  for  m>0  become  unstable  as  ka  increases. 
The  m=l  kink  mode,  on  Figure  1,  becomes  stable  as  ka  nears 
zero.  For  ka*l,  however,  there  is  a  k-band  where  the  m=l 
mode  has  a  larger  growth  rate  than  the  m=0  mode.  The 
linear  theory  predicts  that  all  the  modes  will  become  un¬ 
stable  for  ka>>l.  At  very  short  wavelengths,  however,  the 
instability  is  destroyed  by  small-scale  turbulence  and  mix¬ 
ing  of  the  plasma.  In  practice,  the  largest  growth  rates  are 
expected  for  ka  ~  1. 

A  mere  recent  approach  to  the  numerical  solution 
of  these  problems  has  been  developed  at  the  University  of 
Texas  at  Austin4,  and  consists  of  a  finite  element  technique. 
The  equilibrium  fields  are  developed  in  a  representation 
by  B-splines,  which  form  the  basis  functions  for  the  finite- 
element  solution.  The  differential  equation  forFris  then 
represented  in  difference  form,  for  £  described  as  spline 
coefficients.  The  splines  are  selected  to  automatically 
satisfy  the  boundary  conditions  on  tr-  The  eigenvalue  problem 
is  then  solved  directly,  by  constructing  the  character i Stic 
determinant  ar.d  evaluating  its  root,  _  “ .  This  technique 
has  been  implemented  for  NR!,  on  the  AAV  COR  V  \X  computer 
system . 
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The  code,  EGVPRB ,  which  does  this  problem  is  a 
general  eigenvalue  solver.  It  can  solve  any  eigenvalue 
equation  of  the  form 

A  ( r ,  a  )  f  ~  +  B(r,*K  '  +  C  (r,X)C(r)  »  0 

where  £(r)  is  the  eigenfunction,  X  is  the  eigenvalue  and 
prime  (')  denotes  differentiation  with  respect  to  r.  The 
codes  uses  B-spline5  basis  functions,  which  we  denote  ^.(r). 
Every  function  of  r  is  represented  by  its  spline  fit, 


# 
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A ( t , A )  -  Iiai(X)^i(r) 

B ( r . X )  =  Eibi(X)^i(r) 

C  (r  ,  X)  =  IjC.  (X)U>.(r) 

C(r)  =  (r) 

The  differential  eigenvalue  equation  is  then 
*  *  * 

.a.y.();.ij).  +  2  .  ■  b  •  Y  •  Ui  •  \1>  ■  +  2  -  -C-Y-ib-ib-  =  0. 

ij  i  J  i  J  ij  ij  iyrj 

Multiplying  by  ip  ?  ,  for  each  l  value,  and  averaging 
over  r  (denoted  by  <>)  ,  leads  to  a  matrix  equation, 


VWWf  >  +  +  Y.i  =  °’ 


ft 


or 


3  *  J  3 


or 


rj  •  y  =  o . 

This  equation  has  a  solution  for  y,  provided 

det  M  =  0, 
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which  is  solved  for  the  eigenvalue',  X  =  using  a 

root  -  finder . 

The  code  can  also  be  used  in  a  mode  which  displays 
the  behavior  of  A(r,X),  B(r,A),  and  C(r,A)  vs.  r  as  X  is 
varied,  and  which  shows  det  K  vs.  X.  The  code  is  described 

in  some  detail  in  Appendix  B.  A  simple  test  calulation  is 
described  here  to  illustrate  the  use  of  the  code.  An  equi¬ 
librium  consisting  of  a  uniform  density  plasma  cylinder 
18  -  3 

(n  =  10  cm  ),  of  radius,  r  =  1  cm,  is  assumed  to  carry 
a  uniformly-distributed  current,  I  =  3MA .  The  magnetic 
field  then  rises  linearly  within  the  cylinder,  achieving 
a  peak  value,  Bg(rp)  =  60T  (or  600  KG),  at  the  plasma  edge. 
The  plasma  is  imagined  to  have  a  uniform  temperature,  T  = 
IKeV ,  implying  a  pressure,  p  =  nKgT  =  1.6  *  108  nt/m". 

For  a  kink  mode  with  m  =  1  and  Kr,  =  5,  Figure  3-2 


displays  the  value  of  the  determinant  for  various  trial 

eigenvalues,  where  the  eigenvalues  have  been  normalized  on 

the  plot  so  that  0.194  -  (ojr^/V^)"  S  0.  Here  is  the 

Alfven  speed  at  the  plasma  edge..  The  horizontal  line  is 

det  |M|  =  0,  and  the  intersection  of  the  two  lines  is  the 

root.  Figure  5-5  through  5-5  shows  the  behavior  of  A(u)7  r)  , 
2  2  1 

B(w,  r)  ,  C(u,  r)  as  w  is  varied,  in  ten  equal  steps,  over 
the  range  -  0.194  5  (wr^/V^)^  -  -  0.155,  which  includes  the 
root,  or  to  "manually"  locate  the  root. 

Having  found  an  approximate  root,  the  root  finder 
in  EFVPRB  can  be  utilized  to  refine  the  calculation.  In 
this  case,  the  root  finder  obtained  a  root  at  (wr^/V^)”  = 
-0.185,  in  good  agreement  with  Figure  5-2. 
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Figure  3-3:  Coefficient  A(r 
-0.194<  (wrb/VA 

Kink  Mode:  m= 1 


SECT  I  OX  4 


MAGNETIC  INSULATION 

A  central  issue  in  the  scaling  of  pulsed-power  driv¬ 
ers  to  higher  power  is  whether  the  vacuum  power  feed  to  the 
diode  can  withstand  the  higher  electrical  stress  without 
loss  of  magnetic  insulation.  Magnetic  insulation^  0  refers 
to  the  ability  of  an  applied  magnetic  field  to  turn  emitted 
electrons  back  onto  the  emitting  surface,  thereby  preventing 
electrical  breakdown.  This  concept  is  now  widely-used  in 
the  design  of  high  -power  vacuum  sections  of  pulsed-power 
generators,  vacuum  transmission  lines  and  ion  diodes.^*  The 
required  magnetic  field  can  be  applied  by  external  coils, 
or  it  can  be  the  self-field  due  to  the  current  flowing  across 
the  anode-cathode  gap  at  the  center  of  the  diode. 

Maxwell  Laboratories,  Inc.  has  conducted  a  series  of 
tests  to  study  the  scaling  of  the  loss  current  with  gap  width 
and  electric  and  magnetic  field  strength.  The  apparatus  used 
in  the  experiments  is  shown  schematically  on  Figure  4-1.  It 
consists  of  a  disc  feed  with  a  variable  gap  and  a  5  nH  short- 
circuit  post  with  a  radius  of  5  cm.  The  anode  surface  is 
instrumented  with  a  series  of  Faraday  cup  collectors,  shown 
schematically  on  Figure  4-2,  located  on  a  radial  spoke  at 
various  radii.  The  apertures  of  the  Faraday  cups  are  covered 
with  .003  inch  thick  aluminum  foil  to  shield  the  collector 
from  stray  plasma,  thereby  reducing  the  noise  level. 


47 


r 


!  f  • 


The  experiment  consists  of  driving  current  radially 
into  one  disc  feed  across  the  short-circuit  post,  with 
return  current  flowing  out  on  the  second  disc  feed.  An 
inductive  voltage,  I.dl/'dt,  develops  across  the  gap.  The 
current  flowing  in  the  post  sets  up  an  azimuthal  magnetic 
field  uhich  provides  magnetic  insulation.  The  object  of 
the  experiment  is  to  measure  the  electron  loss  current,  ie. 
the  current  carried  by  free  electrons  which  cross  the 
magnetical ly- insulated  gap.  The  Faraday  cups  provide  a 
measure  of  this  loss. 

The  four  Faraday  cups  used  in  the  experiment  are 
located  at  the  following  locations 


FC#  1 

64cm 

FC#  2 

52cm 

FC#  3 

38cm 

FC#  4 

25cm 

? 

each  having  a  collecting  area  of  1.8  cm  .  The  baffle 
holding  the  foil,  shown  on  Figure  4-2,  limits  the  incident 
angle  for  a  trajectory  to  reach  the  collector.  Ignoring 
scattering  in  the  foil,  only  particles  with  incident  angle, 
9,  relative  to  the  normal  to  the  foil  of  less  than  approxi¬ 
mately  75°  will  reach  the  collector. 


1  f‘ 
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Data  has  been  recorded  for  twenty-one  different 
configurations,  characterized  by  various  gap  widths  and 


A 

driving  currents.  In  this  report  an  analysis  of  one  of 
these  configurations  is  presented.  The  analysis  has  been 
carried  out  with  the  MASK  code,  a  two-dimensional,  fully- 
relat ivist ic  ,  electromagnetic,  particle- in-cel  1  (PIC)  plasma 
simulation  code  developed  by  A.  Drobot  of  SAI  ,  in  collaboration 
with  NRL  ,  MIT  and  Lawrence  Livermore  National  Laboratory. 

The  data  for  this  shot.  Shot  Number  1105,  was 
graciously  supplied  by  John  Shannon  of  Maxwell  Laboratories, 
and  is  summarized  in  Table  4-1.  The  time  history  data  for 
this  shot  is  shown  on  Figures  4-3  through  4-5.  The  drive 
current,  I(t),  is  shown  on  Figure  4-5.  Figure  4-4  is  dl/dt, 

Avhile  Figure  4-5  shows  current,  voltage  and  power. 

Figure  4-6  shows  the  Faraday  cup  and  PIN  diode 
waveforms  for  this  shot.  The  PIN  diodes,  located  adjacent 
to  the  Faraday  cups,  provide  local  x-ray  data,  from  which 
information  about  the  energy  spectrum  of  the  loss  electrons 
can  be  inferred.  In  this. case,  the  PIN  diodes  indicated 
the  presence  of  electrons  with  energy  in  excess  of  100  kev , 
but  the  fraction  of  electrons  above  this  energy  is  not  known 
from  this  diagnostic. 

The  voltage  waveforms  for  the  Faraday  cups  are  the 
voltage  developed  by  the  Faraday  cup  current  as  it  passes 
through  a  50  Q  termination,  with  a  ten-fold  attenuator. 


The  voltage,  Y  ,  indicated  on  the  scope  can  be  translated 
to  Faraday  cup  current,  1^.  ,  by 

10Y  (volts) 

If  (amps)  =  - 

5  on 

The  configuration  shown  in  Figure  4-1,  using  a 
5  cm  gap,  has  been  set-up  with  the  MASK  code,  and  gridded 
on  a  64x16  r-z  mesh.  The  experimental  current  waveform, 
Figure  4-3,  is  used  to  drive  the  simulation.  The  first  test 
of  the  numerical  model  is  a  "cold  test",  which  refers  to 
a  run  in  which  no  free  electrons  are  allowed  to  be  emitted. 
This  type  of  run  tests  the  circuit  model.  Without  current 
smoothing,  the  calculated  induced  voltage  shows  alot  of 
hash.  By  smoothing  the  current  waveform  in  Figure  4-5  to 
relax  sudden  changes  in  the  current,  the  calculated  voltage 
is  found  to  be  in  excellent  agreement  with  the  experimental 
waveform . 

The  cold  test  run  also  provides  a  wealth  of 
information  about  the  field  structure  in  the  device  without 
particles.  Some  of  this  data  is  shown  on  Figures  4-7 
through  4-10,  and  may  be  used  in  conjunction  with  later 
figures  to  examine  the  effect  of  including  the  emitted 
electrons.  On  these  figures,  the  coordinates  are  numbered 
as  X,  =  Z,  X-,=  r,  X,  =  0.  so  that  E,  =  E_  .  E-,  =  E__  and 


m 
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Figure  4-"  shows  the  electric  and  magnetic  field 
spatial  profiles,  as  contours  on  the  r-z  grid,  and  as  a 
vector  plot  for  the  electric  field  lines  (the  magnetic 
field  lines  are  purely  azimuthal).  These  figures  represent 
the  system  at  time  t  =  108  ns,  but  they  are  quite  insensitive 
to  time.  After  the  current  maximum,  dl/dt  switches  sign, 
and  at  later  times  the  electric  field  vector  plot  sho\s-s 
arrows  pointing  opposite  to  those  on  Figure  4-7. 

The  electric  and  magnetic  field  components  Ev  and 
Bfi  are  plotted  against  time  for  various  radii  on  Figure  4-8 
and  4-9.  These  plots  show  the  shape  of  the  driving  current 
waveform  (since  Bg  ®  I)  and  the  shape  of  the  induced  voltage 
waveform  (since  E„  *  V  «  dl/dt). 

Figure  4-10  shows  the  total  field  energy  against 
time.  It  also  displays  the  breakdown  of  field  energy 
associated  with  each  field  component.  The  curves  have 
the  expected  shape  for  an  inductively  driven  gap. 

The  next  step  is  to  turn  on  electron  emission  on 
the  cathode.  To  compare  with  Maxwell's  data,  the  anode 
surface  includes  absorber  regions  which  collect  charge  as 
a  Faraday  cup.  The  absorbers  in  the  code  are  rings  located 
on  the  anode  surface  at  the  same  radii  as  the  Faraday  cups 
in  the  experiment.  Each  absorber  ring  is  four  cells  (or 
4.12  cm)  thick  with  radius  R,  and  therefore  presents  a 
collecting  area  of  2ttP.  x  4.12  cm2,  which  must  be  renormalized 
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to  the  1.8  cm”  collecting  area  of  the  Faraday  cups  used 
in  the  experiments.  The  code  does  not  include  the  .003  inch 
aluminium  foil  or  the  baffle  on  which  the  foil  is  stretched 
(see  Figure  4-2).  All  charges  which  hit  the  absorber  regions 
are  collected  and  counted,  whereas  the  experiment  only 
counted  those  electrons  with  sufficient  energy  to  penetrate 
the  foil  and  with  incident  angles  less  than  approximately 
"5°  to  the  normal. 

The  results  of  the  simulation  are  summarized  on 
Figures  4-11  through  4-24.  Figure  4-11  shows  the  field 
spatial  behavior  at  two  separated  instants  of  time,  t  =  160  ns 
and  t  =  220  ns.  The  earlier  time  shows  a  field  structure, 
which  is  similar  to  the  cold-test  field  (cf.  Figure  4-7), 
except  for  the  effect  of  particles  near  the  short-circuit 
post.  The  late-time  field  structure  shows  the  effect  of 
particle  emission.  Figure  4-12  shows  the  total  field  energy 
and  the  field  energy  associated  with  each  field  component 
plotted  against  time.  These  plots  represent  the  time 
dependence  of  the  volume-averaged  fields.  The  corresponding 
plots  in  the  absence  of  emitted  electrons  are  shown  on 
Figure  4-10.  The  field  energy  is  dominated  by  the  energy 
stored  in  the  Bgf ield . 

The  particle  density  on  the  grid  is  displayed  on 
Figure  4-13  at  t  =  160  ns,  220  ns,  and  340  ns,  which  shows 
the  development  of  the  electron  loss  current  in  the  gap. 


The  total,  volume  -  integrated  charge  in  the  system  plotted 
against  time  is  also  shown  on  Figure  4-13.  The  maximum 
total  charge  in  the  gap  occurs  at  t  =  160  ns,  but  is  highly 
localized  near  the  cathode  surface. 

The  particle  phase  space  projections  are  given  on 
Figures  4-14  through  4-18,  for  time  t  =  160  ns,  240  ns  and 
340  ns.  Figures  4-14  and  4-15  illustrate  the  behavior  of 
the  axial  momentum  vs.  z  (or  X^)  and  vs.  r  (or  X?) , 
and  show  the  development  of  the  electron  loss  current  from 
magnetically- insulated  emission  on  the  cathode  to  the 
formation  of  an  electron  layer  throughout  the  gap,  due 
partially  to  electron  emission  from  the  short-circuit  post. 
Figures  4-16  and  4-17  show  the  same  information  for  the 
radial  momentum,  P^.,  while  Figure  4-18  shows  Pr  vs  P„  . 

The  electron  loss  begins  on  the  cathode  at  large  r,  where 
the  magnetic  insulation  is  weakest,  and  gradually  progresses 
toward  the  short-circuit  post,  since  the  current  crossing 
the  gap  at  large  r  reduces  the  magnetic  field  at  small  r. 

The  E'J  instantaneous  power  is  plotted  against  time 
on  Figure  4-19,  which  shows  the  total  power,  as  well  as 
the  contributions  due  to  the  radial  and  axial  components. 

The  axial  component,  E,JZ,  is  the  dominant  one,  as  expected 
for  current  loss  across  the  gap. 

Electron  emission  in  the  code  is  allowed  from  both 
the  cathode  surface  (or  right-hand  boundary)  and  the  surface 


of  the  short-circuit  rod  (or  lower  boundary).  Figure  4-20 
shows  the  emitted  current  from  these  two  surfaces  plotted 
against  time. 


I 


I 


r> 


# 


!  f* 


The  electrons  may  be  collected  on  all  surfaces. 
Figures  4-21  shows  the  collected  current  vs.  time  on  the 
anode  (left),  cathode  (right),  top  (upper)  and  short-circuit 
rod  (lower)  surfaces. 

The  instantaneous  current  collected  by  absorber 
number  1  through  4 ,  which  correspond  to  Faraday  cups  1 
through  4,  is  shown  on  Figure  4-22.  The  total  integrated 
current,  and  the  current  integrated  over  1000  time  steps 
(10  sec)  is  shown  on  Figure  4-23  for  Faraday  cup  number  1 
and  on  Figure  4-24  for  Faraday  cup  number  2. 

The  calculated  current  on  Figure  4-24  has  two  major 
discrepancies  with  the  experimental  data.  First,  the 
waveform  for  the  current  averaged  over  10  ns  bins  shows 
two  peaks  separated  by  almost  200  ns  in  time.  The  experi¬ 
mental  Faraday  cup  waveform  does  not  show  such  widely 
separated  peaks.  Second,  the  magnitude  of  the  integrated 
loss  as  calculated  by  the  code  is  approximately  600  times 
greater  than  the  experimentally-measured  loss. 

Both  of  these  discrepancies  can  be  attributed  to 
the  absence  in  the  code  of  the  .003  inch  aluminium  foil 
covering  the  Faraday  cup  and  the  geometrical  aperture 
caused  by  recessing  the  Faraday  cup  charge  collector  as 
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shown  on  Figure  4-2.  The  aluminium  foil  will  stop  a  130  keV 
electron  at  normal  incidence.  A  grazing  electron  with  much 
lower  energy  will  stop  in  the  foil. 

An  examination  of  the  particle  energy  spectrum 
computed  by  the  code  indicates  that  essentially  all  of  the 
electrons  which  make  up  the  first  spike  on  Figure  4-24  lie 
below  100  keV  and  therefore  would  not  be  detected  in  the 
experiment.  Half  of  the  electrons  in  the  second,  later 
spike  on  Figure  4-24  also  lie  below  100  keV.  The  other 
half  are  energetic  enough  to  produce  the  PIN  diode  signal 
observed  in  the  experiment. 

The  electron  orbits  as  they  impinge  on  the  absorbers 
(Faraday  cups)  in  the  code  are  very  steep  since  the  electrons 
drift  radially- inward  as  they  traverse  the  gap,  due  to  the 
Ez  x  B0  drift.  Most  of  these  electrons  are  therefore 
blocked  by  the  Faraday  cup  acceptance  (0<  75°)  aperture, 
or  are  stopped  in  the  aluminium  foil. 

Modifications  to  MASK  are  currently  in  progress  to 
quantify  these  effects.  The  results  to  date,  however, 
indicate  that  the  electron  losses  measured  in  the  Maxwell 
experiment  probably  constitute  only  a  small  fraction  of 
the  actual  loss  present  in  the  apparatus.  The  scaling  of 
existing  pulse-power  devices  to  significantly  higher  power 
will  depend  on  the  understanding  and  control  of  these  power 
losses  in  the  vacuum  sections  of  the  machines. 
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MAXWELL  DATA  FOR  SHOT  1105 


T  <01  > 


Experiment;)  1  dl/dt  Waveform 


Measured  Current  (Squares),  Voltage  (Triangles) 
and  Power  (Solid)  Waveforms 


MAXWELL  SHOT  NO.  1106 


FIGURE  4-6:  FARADAY  CUP  AND  PIN  DIODE  WAVEFORMS 
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Fig.  4-8:  Cold  Test  Ez  v£.  time  on  the  Anode  Surface 

(a)  r=  21.5  cm;  (b)  r=38.0  cm; 

(c)  r«  52.4  cm;  (d)  r=63.7  cm. 


Fig.  4-10:  Cold  Test  Field  Energy  vs^  Time. 

(a)  Dz  Component;  (b)  D  Component; 

(c)  HQ  Component;  (d)  Tital  Field  Energy 
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Fig.  4-12:  Stored  Field  Energy  vs.  time. 

(a)  D  Component;  (b)  D_  Component; 

2  r 

(c)  Hq  Component;  (d)  Total  Field  Energy. 
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Fig.  4-13:  Particle  Density  (r-z  plots)  and  Total  Charge  on  the  Grid  vs. 
time. 


(a)  Particle  Density  at  t*160  ns 

(b)  Particle  Density  at  t*220  ns 

(c)  Particle  Density  at  t-340  ns 

(d)  Total  Charge  on  Grid. 
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4-14:  Axial  Momentum,  P„ ,  vs.  Z. 

z  — 

(a)  t*160ns;  (b)  t=220ns;  (c)  t*340ns, 
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4-16:  Radial  J 
(a)  t*160ns;  (I 
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Fig.  4-21:  Collected  current  vs.  time. 

(a)  Current  collected  on  left  surface; 

(b)  Current  collected  on  right  surface; 

(c)  Current  collected  on  lower  surface; 

(d)  Current  collected  on  upper  surface. 


Appendix  A 
WIRES  CODE 


The  WIRES  code,  based  on  the  model  described  in  Section 
2,  is  basically  a  Runge-Kutta  integrator  for  five  variables: 

(1)  array  radius, 

(2)  implosion  speed, 

(3J  current, 

(4)  radiation  yield  in  photons  with  energy  greater 
than  e*  (an  input), 

(5)  total  radiation  yield. 

If  run  interactively,  the  code  will  prompt  the  user  for 
the  following  data: 

Block  1 

N  =  number  of  wires 

EST=  spectrum  cut-off  energy  (eV) 

XMU=  single  wire  mass  per  unit  length  (g/cm) 

XL  =  wire  length  (cm) 

Block  2 

R  ( 0 )  - 

B 
Z 

XMASS  = 

CLOG  = 

GAMMA  = 

EMISS1= 

EMI SS2= 

NPFLAG= 


Block  5 

V0  =  circuit  charge  voltage  (Volts) 
Z0  =  generator  impedance  (Ohms) 

XLD=  diode  inductance  (Henries) 


initial  array  radius  (cm) 

outer  radius  for  return  current 

atomic  number  of  wire  material 

atomic  mass  (amu)  of  wire  material 

Coulomb  logarithm  (default  value  =  4) 

specific  heat  ratio  (default  value  =  5/3)  _4 

emissivity  for  hv<EST  (default  value  =  5.x  10,) 

emissibity  for  hv>EST  (default  value  =  5.x  10  ) 

(1  or  0)  =  (Yes  or  No)  print  during  integration 

for  wire  cooling  after  assembly. 


DT  =  time  step  for  Runge-Kutta  (sec) 

NPRINT  =  number  of  time  steps  between  print-out's. 


Each  data  block  should  be  entered  in  free-format  as  a  single- 
line  input. 


The  main  program  initializes  the  problem  and  calls  the 
following  subroutines: 


(1)  STEP:  Calculates  one  Runge-Kutta  time  step, 

using  subroutine  FORCE  to  calculate  the 
necessary  first  derivatives. 


(2)  FORCE:  Provides  derivatives  for  use  by  STEP. 

FORCE  finds  the  temperature  by  imposing 
a  local  Bennett  equibrium, 

2 

5?  '  "  (l*Zeff)KBT. 

(3)  RADFRAC:  Calculates  the  fraction  of  the  black- 

body  radiation  yield  which  lies  above 
EST. 


(4)  XCURR:  Allows  a  specified  current  waveform  to 

be  utilized;  this  option  is  not  used  in 
the  current  version  of  WIRES. 


(5)  OUT:  Print-out  subroutine.  The  following 

quantities  are  printed: 


(a) 

T 

= 

time 

(b) 

Y(l) 

= 

array  radius 

(c) 

Y  (2) 

= 

implosion  speed 

(d) 

Y(3) 

= 

current 

(e) 

Y  (4) 

yield  for  hv>EST 

(f) 

Y  (5) 

= 

total  yield 

(g) 

A 

= 

wire  radius 

(h) 

T  (EV) 

= 

wire  temperature  (eV) 

(i) 

ZEFF 

= 

effective  ionization  state 

O) 

DENS 

= 

number  density 

00 

RP (OHMS) 

- 

wire  resistance (Spitzer) 

(1) 

LP(H) 

= 

wire  inductance  (Russell) 

(m) 

VI  (W) 

= 

ipput  power  =  IV 

(n) 

P-OHMIC 

= 

I  Rp  =  Ohmic  dissipation 

(o) 

P-BB 

= 

blackbody  radiation  power 
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(p)  K-FLD  =  LI  / 2  =  stored  field  energy 

(q)  W-KIN  =  N p £ v  /2  =  kinetic  energy 

(rj  W-INT  =  1.5  n(l+Zg££)  KgT  =  internal  energy 

(6)  XVOLTS:  specifies  applied  voltage  waveform. 

(7)  FINAL:  Calculates  final  assembly  and  cooling  of 

plasma  cylinder,  based  on  instantaneous 
conversion  of  kinetic  energy  to  temperature 
followed  by  radiative  cooling  via  blackbody 
emission. 

(8)  DERIV:  Provides  derivatives  for 

(a)  temperature 

(b)  yield  above  EST 

(c)  total  radiation  yield 

needed  by  subroutine  FINAL. 


Appendix  B 
EGVPRB  CODE 

The  EGVPRB  Code,  together  with  pre -processors  (EGVSETUP 
and  MHDEQUIL)  and  post-processors  (EGVPLT  and  EGCOPLT)  are 
described  in  an  on-line  documentation  file,  EGVPRB. INF,  which 
is  listed  below. 

The  various  modules  are 

(1)  EGVSETUP  :  File  assignment. 

(2)  MHDEQUIL. FOR  :  MHD  equilibrium  specification. 

(3)  EGVPRB. FOR  :  Linear,  ideal  MHD  stability  analysis. 

(4)  EGVPLT. FOR  :  Plots  coefficients  and  eigenfunction 

for  converged  solution. 

(5)  EGCOPLT • FOR  :  Plots  coefficients  and  eigenfunctions 

for  sequence  of  trial  solutions. 

The  main  module,  EGVPRB. FOR,  solves  general  second-order, 
differential  eigenvalue  problems  of  the  form, 

a +  bC"  +ce  •=  0, 

where  £(r)  is  'the  eigenfunction,  and  the  coefficients,  a(r,X), 
b(r,X),  and  c(r,X),  are  functions  of  both  r  and  the  eigenvalue 
parameter,  X.  This  code ,  speciali zed  to  solve  the  linear,  idea 
MHD  stability  problem  for  a  specified  cylindrical  equilibrium, 
is  set-up  on  the  JAYCOR  VAX  Computer  (Host  CAIN,  Directory 
£lPR3j).  The  MHD  stability  analysis  itself  is  described  in 
Section  3. 


The  maj or  subunits  o f  the  EGVPRB  code  are  the  following: 


(1) 

SYSODE : 

the  main  program,  a  top-level  governor. 

(2) 

MATRIX: 

reads  input  data  and  calculates  the  a,b,c 
coe  f  f icients . 

(3) 

FCN : 

calculates  the  determinant  which  provides 
the  characteristic  equation  for  the  eigen 
value . 

(4) 

CEVALF: 

a  root  finder  which  solves  the  character¬ 
istic  equation  for  the  eigenvalue. 

(5) 

BC: 

a  subroutine  which  sets-up  the  specified 
boundary  conditions. 

(6) 

MATNRM : 

normalizes  the  determinant  to  avoid  over¬ 
flow/underflow. 

(7) 

NRMFCN : 

eigenfunction  calculation 

(8) 

PLTFCN : 

sets-up  output  plots 

(9) 

DEPSE : 

decomposition  of  function  into  B-splines. 

(10) 

REPSE: 

recomposition  of  function  from  B-splines. 

(ID 

REPSP: 

recomposition  of  first  derivative  of 
function  from  B-splines. 

The  EGVPRB  package  was  designed  for  the  CRAY  computer 
system,  and  several  CRAY-dependent  lines  of  code  were  "com- 
mented-out"  of  the  code  to  adapt  it  to  the  VAX.  The  unfor¬ 
tunate  overflow/underflow  limits  on  the  VAX  (approximately 
±35 

10  )  impose  a  limitation  on  the  number  of  knots  (or  nodes) 

which  may  be  carried  in  the  splines.  The  determinant  to 
be  computed  is  an  N  *  N  determinant,  where  N  is  the  number 
of  knots.  While  the  determinant  is  normalized,  the  VAX  can 
overflow  or  underflow  easily  if  N  t  30  is  utilized.  An 


input  parameter,  SETNRM,  has  been  built  into  the  code  to 
"fine  tune"  the  determinant  normalization  so  as  to  avoid 
this  problem.  With  SETNRM  specified  as  1,  the  code  normal 
izes  the  determinant  to  the  largest  element  in  the  matrix. 
Test  problems  using  N=20  have  run  without  difficultly  with 
SETNRM=1  . 

Listings  of  the  various  modules  follow. 


APPENDIX  C 


sting  of  WIRES 
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PROGRAM  WIRES 
IMPLICIT  REAL*8 (A-H,0-Z> 

PARAMETER (ND I M=5) 

DIMENSION  Y(NDIM) ,DY(4,NDIM) >  YOLD (NDIM) 

COMMON  N,EST,  XMU,RH0,  XL, B, Z, CLOG, NPFLAG, 

*  X MASS ,  GAMMA , DT , TEMP , A , EM I SS 1 , EM I SS2 
COMMON/CIRC/VO, ZO, XLD, XLDOT , XLP.RP, ZEFF 
DATA  PI/3. 141 59265358979323BD0/ 

PROGRAM  TO  CALCULATE  IMPLOSION 
OF  WIRE  ARRAYS 
Y(1)=ARRAY  RADIUS 
Y  <  2 ) — I MPLOS I ON  SPEED 
Y (3) =CURRENT 
Y (4) -YIELD  ABOVE  EST 
Y<5)=T0TAL  YIELD 

INPUTS 

N=NUMBER  OF  WIRES 
EST=CUT-OFF  ENERGY  (EV) 

XMU=WIRE  MASS/LENGTH 
XL=WIRE  LENGTH  (CM) 

R (0) = INITIAL  ARRAY  RADIUS  (CM)  =Y ( 1 ) 

B=OUTER  RADIUS  FOR  RETURN  CURRENT  (CM) 
Z=ATOMIC  NUMBER 
CLOG=COULOMB  LOG  <DEF:4) 

XMASS=ATOMIC  MASS  (AMU) 

G AMMA=SPEC I F I C  HEAT  RATIO  (DEF:5/3) 

EM I SS 1 =EM I SS I VI T Y  BELOW  EST  (DEFs  5.E-4) 

EM I SS2=EM I SS I V I TY  ABOVE  EST  (DEF:  5.E-6) 

NPFLAG=(1,0)=(YES,N0)  PRINT  DURING  TEMP  DECAY 

VO=CIRCUIT  CHARGE  VOLTAGE 

ZO=GENERATOR  IMPEDANCE 

XLD=DIODE  INDUCTANCE 

DT=TIME  STEP  (SEC) 

NPR I NT= INTERVAL  BETWEEN  PRINTS 


GAMMA=5. /3. 

CL0G=4. 

EM I SSI =5. E-4 
EMISS2=5. E-6 
Y (2) =0. 

Y (3) =1 . E5 
Y (4) =0. 

Y (5) =0. 

IFIRST=0 
PRINT  900 

FORMAT (IX, ' N, EST (EV) , XMU, XL' / 

*  IX, 'R(O) , B , Z , XMASS ( AMU ) , CLOG, GAMMA, ' , 

* ’EMISS1,EMISS2, NPFLAG ' / 

*  5X,'NPFLAG=1  FOR  PRINT'/ 

*  IX, ' VO(VOLTS) , ZO (OHMS) , XLD (HENRIES) ' / 

*  IX, ' DT , NPR I NT' ) 

READ  <5, # )  N, EST , XMU, XL 

READ (5,*)  Y ( 1 ) , B, Z , X MASS , CLOG , GAMMA , EMI SSI, EM I SS2 , NPFLAG 
READ (5,*)  VO, ZO, XLD 
READ  <5, #)  DT, NPR I NT 
T=0. 

K0UNT=0 

KPRINT=0 

CALL  STEP(T, Y, NDIM, DY, YOLD) 

K0UNT=K0UNT+1 


K  PR I NT=K  PR I NT+ 1 
XN=N 

IF  < ( I FIRST . EQ. 0) . AND.  (A.GT.Y(l) #SIN (PI /XN) ) )  GOTO  1 
IFIRST=1 

IF (ABS(A-Y (1) *SIN(PI/XN) ) .LT. l.E-3*A)  GOTO  20 
DT=MIN(DT, . 5* (A/SIN (PI /XN) -Y ( 1 ) ) /Y (2) ) 

IF (KPRINT.LT. NPRI NT)  GOTO  10 
KF’RINT=0 

CALL  OUT  <T,Y,NDIM) 

GOTO  10 

CALL  OUT (TiYt ND I M ) 

CALL  FINAL <T,Y, NDIM) 

STOP 

END 

SUBROUT I NE  STEP ( T , Y , ND I M , DY , YOLD ) 

IMPLICIT  REALMS ( A— H,  0— Z ) 

DIMENSION  Y(NDIM) , DY(4,NDIM) , YOLD (NDIM> ,D(4> 

COMMON  N,EST,  XMU,RHO,  XL ,  B,  Z ,  CLOG,  NF'FLAG, 

*  XMASS , GAMMA , DT , TEMP , A , EM I SS 1 , EM I SS2 
TOLD=T 

DO  5  1=1, NDIM 
YOLD ( I ) =Y ( I ) 

D(l)=DT/2. 

D (2) =DT/2. 

D(3)=DT 

D(4)=DT/6. 

L=1 

CALL  FORCE (T,  Y,  DY, NDIM, L> 

L=L+1 

IF(L.EQ.S)  GOTO  20 
T=T0LD+D(L-1) 

DO  15  J=1 ,  NDIM 

Y ( J) =YOLD ( J) +DY (L-l , J) *D(L-1> 

GOTO  10 

DO  25  J=1 , NDIM 

Y ( J) =YOLD ( J) +D (4) *  <DY(1,  J>+2.  *DY(2,  J) 

*  +2.  #DY (3, J)  +DY  (4,  J)  ) 

RETURN 

END 

SUBROUT I NE  FORCE ( T , Y , DY , ND I M , LRK > 

IMPLICIT  REAL*8(A-H,0-Z> 

DIMENSION  Y (NDIM) »  DY (4*  NDIM) 

COMMON  N , EST , XMU , RHO , XL , B , Z , CLOG , NPFL AG , 

*  XMASS , GAMMA , DT , TEMP , A , EM I SS 1 , EM I SS2 
COMMON/CIRC/VO, ZO, XLD, XLDOT, XLP, RP, ZEFF 
DATA  P I / 3 . 1 4 1 592653589793238D0 / 

DATA  SIG/5.6696E-5/, XKB/1 . 3B07E-16/ 

DATA  TO/1. 16E7/, XMP/1.6606E-24/ 

Z26= (26. /Z) **2 

XN=N 

XIP=Y (3) 

C=XIP*XIP#XMASS*XMP/200.  /XN**2/XMU/XKB 
TU=C 

TL=C/ (l.+Z) 

DO  10  1=1,20 
TEMP=. 5* (TU+TL) 

ZEFF=26. *SQRT (TEMP/ (T0+Z26CTEMP) ) 

CTEST=TEMP# ( ZEFF+ 1 .  ) 

IF (CTEST. LT. C)  GOTO  5 

TU=TEMP 

GOTO  10 

TL=TEMP 

CONTINUE 


IF ( ZEFF . LT. 2)  GE=.582+. 101  * ( ZEFF-1 ) 

IF ( ( ZEFF. LT. 4) .AND.  (ZEFF.GE.2) )  GE=. 683+. 051  * ( ZEFF-2) 

IF ( ( ZEFF . LT. 16) .AND. (ZEFF. GE. 4) ) GE=. 785+. 01 15# ( ZEFF-4) 

IF (ZEFF. GE. 16)  GE=1 . -1 . 232/ZEFF 
RH0=3800. #ZEFF*CLDG/GE/TEMP## ( 1.5) 

CALL  RADFRAC  (EST »  TEMP,  FRAC) 

B:!MISS=EMISS1  #  ( 1 .  -FRAC)  +EMISS2*FRAC 

A=( l.E7*RH0#XIP*XIP/2. /PI/PI /XN/XN/SIG/EMISS/TEMP*#4) **( 1 . /3. ) 
RP=RH0*XL/PI/A**2/XN 

XLP=XL# ( . 5+2. #L0G (B##N/XN/A/Y ( 1 ) ## (N-l ) ) ) /XN#1 . E-9 
XLD0T=— 2. #XL# (XN-1 . ) #Y (2) /Y ( 1 ) /XN#l.E-9 
CALL  XV0LTS(T,V) 

AS=2. #PI #A#XL#XN 
DY (LRK, 1 ) =Y (2) 

DY (LRK, 2) =- (XN-1. ) * (XIP/10. /XN) ##2/XMU/Y ( 1 ) 

DY (LRK , 3) = (V- ( Z0+RP+XLD0T) *X IP) / ( XLD+XLP) 

DY(LRK,4)=FRAC#SIG*EMISS2#AS*TEMP**4 

DY (LRK , 5) =SIG*EMISS#AS*TEMP**4 

RETURN 

END 

SUBROUT I NE  RADFRAC ( EST , TEMP , FRAC ) 

IMPLICIT  REAL #8 ( A— H, 0— Z ) 

DIMENSION  Y'LT  (70)  ,  YFRAC  (70) 


DATA 

YLT/ 

.01, 

.02, 

.03, 

.04, 

.05, 

.055 

* 

.065 

i,  .07 

,  .07 

5  7 . 03 7  ■  085 7  a 

09,  . 

095, 

% 

.  10, 

.11, 

.  12, 

.  13, 

.  14, 

.  15, 

.  16, 

.  17, 

% 

.20, 

.22, 

.24, 

.26, 

.28, 

.30, 

nn 

•  Ui.  7 

.34, 

* 

.40, 

.45, 

.50, 

.  55, 

.60, 

.65, 

.7,  . 

8,  .9 

% 

1.1, 

1.2, 

1.3, 

1.4, 

1.5, 

1.6, 

1.7, 

1.8, 

% 

2.5* 

3.  ,3 

-.5,4 

■  7  5i 

7  <S  m  7 

7.  ,8 

.  ,9. 

,  10. 

* 

20.  7 

30.  , 

40.  , 

wjO  a  7 

100. 

/ 

DATA  YFRAC/O. , 3. 7E-27 , 2. 7E-17, 1.9E-12, 
t  1.3E-9, 1 . 35E-8 , 9 . 29E-8 , 4 . 67E-7 , 1 . 84E-6 , 

K  5.94E-6, 1.64E— 5,3.99E-5,8.7E-5, 1.73E-4, 

H  3.21E-4,9. 11E-4, .00213, .00432, .00779, 
t  .01285, .01971, .02853, .03933, .05210, 

*  .06672, . 10087, . 14024, . 18310,-22787, 

*  .27320, .31807, .36170, .40327, .44334, 

*  .  48084 ,  .  56428 ,  .  63370 , ..  69086 ,  .  73777 , 

*  .77630, .80806, .85624, .88998,  .91415, 

*  .93184,. 94505 , . 95509 , . 96285 ,  .  96893 , 

*  .97376, .97765, .98081, . 98340, . 98555 , 

*  .99216, .99529, .99695, .99792,-99890, 

*  . 99935 , . 99959 , . 99972 , . 99980 , . 99985 , 

*  . 999955, . 99998, . 9999943, . 9999975, . 9999988, 

*  .99999985/ 

DATA  HC/ 1 . 2399E-4/ , C2/ 1 . 43BB3/ 

XLT=TEMP*HC/EST 
IF(XLT.GE. .02)  GOTO  10 
FRAC=0. 

RETURN 

IF(XLT.LE. 100. )  GOTO  20 
X=C2/XLT 

FRAC=1.-.0513*X#*3 

RETURN 

DO  30  1=2,70 

IF (XLT. GT. YLT ( I > )  GOTO  30 

FRAC=YFRAC ( I - 1 ) + ( YFRAC ( I ) -YFRAC ( I - 1 ) ) * 

*  ( XLT-YLT ( I— 1 ) ) / (YLT(I)-YLT(I-l) ) 

GOTO  40 

CONTINUE 

RETURN 

END 

'vlvlvlvlv/ v' \ •. 1 


SUBROUTINE  XCURR(XIP,T) 

IMPLICIT  REAL*8(A-H,0-Z> 

X IP=5. E6 

RETURN 

END 

SUBROUTINE  OUT (T>V> NDIM) 

IMPLICIT  REAL*8(A-H,0-Z> 

DIMENSION  Y(NDIM) 

COMMON  N , EST , XMU , RHO , XL , B , Z , CLOG , NPFLAG , 

*  XMASS, GAMMA, DT, TEMP, A, EMISS1 ,EMISS2 
COMMON/CIRC/VO, ZO, XLD, XLDOT, XLP,RP, ZEFF 

DATA  PI/3. 14 159265358 9793238 DO/ , XMP/1 . 6606E-24/ 
DATA  SIG/5. 6696E-5/ , XKB/1 . 38C7E-16/ 

WRITE (6, 999) 

999  FORMAT  (19X,  ’T',11X,  'Y(l)  ',  11X,  'Y(2)  MIX,  *Y(3)  ', 

*  11X, ' Y(4) ' , 11X, ' Y<5) ' ) 

998  FORMAT (34X, 'A' , 10X, 'T(EV) ' ,  11 X ZEFF ', 1 ■ X ,' DENS ’ > 

997  FORMAT (27X, 'RP(OHMS) ' , 10X, 'LP(H) ' , 10X, 'VI (W) ' ) 

996  FORMAT (28X,  'P-OHMIC'  ,  1 1  X  ,  'P-BB'  ,  10X,  'F'-KIN'  ) 

995  FORMAT (30X, ' W-FLD ' , 10X, 'W-KIN' , 10X , ' W-INT ' ) 

TTTT=TEMP/1 1600. 

WRITE (6, 1000)  T, (Y(I) , 1=1, NDIM) 

WRITE (6, 998) 

DENS=XMU/PI/A**2/XMASS/XMP 
CALL  XVOLTS (T, V) 

CALL  RADFRAC (EST, TEMP, FRAC) 

EMISS=EMISS1* ( 1 . -FRAC) +EMISS2*FRAC 
VI=V*Y(3> 

XI2R=RP*Y(3) **2 
XN=N 

AS=2.*PI*A*XL*XN 
PBB=SIG*EMISS*TEMP**4*AS*1 . E-7 
XI2LD=.5*XLD0T*Y(3) **2 
WFLD=.5* (XLD+XLP) *Y(3) **2*1.E7 
WKIN=.5*XMU*XL*Y(2) **2*XN 
WINT=1.5*DENS*  (ZEFF+l'.  )  *XKB*TEMP 
WRITE (6, 1001)  A, TTTT , ZEFF , DENS 
WRITE (6, 997) 

WRITE (6, 1002)  RP, XLP, VI 
WRITE (6, 996) 

WRITE(6, 1002)  X 1 2R , PBB , X 1 2LD 
WRITE (6, 995) 

WRITE(6, 1003)  WFLD , WK I N , W I NT 

1000  FORMAT (5X, 6E15. 5) 

1001  FORMAT (20X,4E15.5) 

1002  FORMAT (20X.3E15. 5) 

1003  FORMAT (20X , 3E15. 5/ ) 

RETURN 

END 

SUBROUTINE  XVOLTS (T,V) 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/C IRC /VO, ZO, XLD, XLDOT, XLFsRP, ZEFF 

v=vo 

RETURN 

END 

SUBROUTINE  FINAL (T, Y, NDIM) 

IMPLICIT  REAL f 8 ( A-H, 0-Z ) 

DIMENSION  Y(NDIM) , ZVECT (3) , DZV (4, 3) , D (4) , Z0LD(3) 
COMMON  N, EST , XMU, RHO, XL, B, Z , CLOG, NPFLAG, 

*  XMASS , GAMMA , DT , TEMP , A , EM I SS 1 , EM I SS2 
COMMON/CIRC/VO, ZO, XLD, XLDOT, XLP.RP, ZEFF 
DATA  P I / 3 . 1 41592653589793238D0/ 

DATA  XMP/1 . 6606E-24/ , TO/ 1 . 16E7 / 
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DATA  S1G/5.6696E-5/,  XKB/ 1 . 3807E-16/ 

FINAL  ASSEMBLY  —  CONVERT  TO  CYLINDER 
VA=ALFVEN  SPEED 

XKK=l/RO=WAVENUMBER  FOR  MAXIMUM  GROWTH 
TASBLY=5#MHD  GROWTH  TIME 
KINETIC  ENERGY  CONVERTED  TO  TEMPERATURE 
RADIATION  ASSUMED  BLACK. -BODY 
ZVECT  ( 1 )  =TEMF'ERATURE 
ZVECT  <2)=RADIATED  ENERGY  ABOVE  EST 
ZVECT (3)=T0TAL  RADIATED  ENERGY 
XN=N 

126- (26. /I) **2 
ROUT=A+Y ( 1 ) 

DENS=XN*XMU/PI/R0UT*#2 

BOUT=Y (3) /5. /ROUT 

VA=SQRT (BOUT  * BOUT/4. /PI /DENS) 

XKK=1. /ROUT 
TASBLY=5. /XKK/VA 
WKIN=.5*XN*XMU*XL*Y(2) **2 
TIN=TEMP 

C=XMASS*XMP#Y (2) **2/ (3. *XKB) 

TU=TIN+C 
TL=TIN+C/ (l.+Z) 

DO  10  1=1,20 
TEMP=.5* (TL+TU) 

ZEFF=26. *SQRT (TEMP/ (T0+Z26*TEMP) ) 

CTEST= (TEMP-TIN) *(1.+ZEFF) 

IF(CTEST.LT.C)  GOTO  5 

TU=TEMP 

GOTO  10 

TL=TEMP 

CONTINUE 

WRITE (6, 998)  TEMP/ 1 1600. , ZEFF 
FORMAT (10X,  ' TEMP, ZEFF=  * , 2E15. 5/) 

CALL  RADFRAC ( EST , TEMP , FRAC ) 

EM I SS=EM I SS 1  * ( 1 . -FRAC ) +EM I SS2  *FR AC 
A5=2. *PI*ROUT*XL 
ZVECT ( 1 ) =TEMP 
ZVECT (2) =Y (4) 

ZVECT (3) =Y (5) 

XK2=SIG*AS 

XK1  =  XK2*XMASS*XMP/  (1.5»XNH:XMU*XL*XKB) 

DT=.05* (l.+ZEFF) /XK1/TEMP**3/EMISS 
D(l)=DT/2. 

D(2)=DT/2. 

D (3) =DT 
D(4)=DT/6. 

KPRINT=0 

NMAX=1+INT (TASBLY/DT) 

NPR I NT =NMAX / 1 00+ 1 
DO  40  JJJ=1 , NMAX 
TOLD=T 

DO  15  KK.K=1 , 3 
ZOLD (KKK) =ZVECT (KKK ) 

L=  1 

CALL  DER I V ( T , ZVECT, DZV,L, XK 1 , XK2, Z26, EST, EMISS1 , EMISS2) 
L=L+1 

IF(L.EQ.S)  GOTO  30 
T=TOLD+D (L-l ) 

DO  25  KKK-1.3 

ZVECT  (KKK. )  =ZOLD  (KKK )  +DZV„(L- 1 ,  KKK )  #D  <  L- 1 ) 

GOTO  20 
DO  35  KKK=1 , 3 


V.V.V.N 


•.U>.  -- 


ZVECT (KKK)=ZOLD (KKK >  +D (4) * (DZV < 1 , KKK) +2. *DZV (2, KKK ) 

*  +2. *DZV(3,KKK)+DZV(4,KKK) ) 

IF (ZVECT ( 1 ) . LE. 1 . E3)  GOTO  50 
I F  ( NF'FL AG .  EG .  0 )  GOTO  40 
K PR  I  NT=K  F’R  I  NT  + 1 
IF (KPRINT.LT. NPRINT)  GOTO  40 
KPR I NT=0 

WRITE (6, 999)  T, ZVECT ( 1 ) / 1 1600. , ZVECT (2) *1 . E-7, ZVECT (3) *1 . E-7 

999  FORMAT ( 1 OX , * T , TEMP , WRADG , WRAD= * , 4E 1 5 . 5 ) 

40  CONTINUE 

50  CONTINUE 

m  TTTT=ZVECT(1) / 11600. 

WRADG* ZVECT < 2) *1. E-7 
WRAD= ZVECT (3) Cl .E-7 

ZEFF*26. CSORT (ZVECT <  1 )  / (T0+Z26*  ZVECT ( 1 ) ) ) 

WRITE (6, 1000)  ROUT, TASBLY, DENS, TTTT, ZEFF , WRADG, WRAD 

1000  FORMAT (1H0,20X, 'FINAL  ASSEMBLY'/ 

m  #  5X, 'COLLAPSE  RADIUS (CM) = ' , E15. 5/ 

•  *  5X, 'ASSEMBLY  TIME (SEC) = ' , El5. 5/ 

*  5X , ' DENSITY (G/CC)  = ' , E15.  5/ 

*  5X, 'TEMPERATURE (EV)*' ,E15. 5/ 

*  5X, ' ZEFF* ' , E15. 5/ 

*  5X, 'RADIATION  ABOVE  EST  <J)=',E15.5/ 

*  5X, ' TOTAL  RAD I AT I ON (J)=' ,E15.5) 

RETURN 

END 

SUBROUTINE  DER I V ( T , ZVECT, DZV,L, XK1, XK2, Z26,EST, 

*  EM I SS 1 , EM I SS2 ) 

IMPLICIT  REALC8 (A— H, 0— Z) 

DIMENSION  ZVECT (3) »DZV(4,3) 

•  DATA  T0/1.16E7/ 

IF ( ZVECT ( 1 ) . LT. 0. )  ZVECT ( 1 ) *0. 

TEMP* ZVECT ( 1 ) 

TFAC=T0+Z26*TEMP 

CALL  RADFRAC (EST, TEMP, FRAC) 

EMISS=EMISS1* ( 1 . -FRAC) +EMISS2CFRAC 

•  DZV(L, 1>=-EMISS*XK1*TEMP**4/ (1.+26. CSGRT (T/TFAC) * ( 1 . +. 5CT0/TFAC) ) 
DZV (L, 2) *EMISS2*FRAC*XK2*TEMP**4 
DZV (L, 3) =EMISS*XK2*TEMP**4 
RETURN 
END 

« 

Ck 


1% 


<« 


EGVPRB.  INF 


This  is  a  user -info  file  for  usiris  the  EGVF'RB  package 
for  solving  eigenvalue  problems  of  the  form, 

ay''  +  by'  +  cy  =  0  , 

where  the  coefficients,  a,b,c,  depend  on  the  independent 
variable,  r,  and  on  the  eigenvalue.  The  code  finds  the 
eigenvalue  as  well  as  the  eigenfunction,  v(r).  The  current 
version  (as  of  8/1/83)  of  EGVF’RB  is  set  up  to  solve  the  linear, 
ideal  MHD  problem  for  an  arbitrary  cylindrical  equilibrium. 

The  user  mav  use  this  package  to  solve  other  eigenvalue  problems 
of  the  form  given  above  bv  defining  new  coefficients  and 
boundary  conditions  in  SUBROUTINE  MATRIX. 

The  code  may  be  run  in  either  of  two  modes.  In  the  first 
mode  it  calls  COMPLEX  FUNCTION  CEVALF  to  find  the  roots 
of  thh  characteristic  eigenvalue  equation.  Alternatively, 
the  code  may  be  used  in  a  mode  where  it  examines  the  value 
of  the  eigenvalue  equation  over  a  range  of  user-specified 
trial  eigenvalues.  This  second  mode  allows  the  user  to  search 
manually  for  the  root,  and  to  examine  the  behavior  of  the 
coefficients  as  the  eigenvalue  parameter  is  varied. 


STEP  1 

©esvsetup 

This  command  causes  the  VAX  to  assign  names  to  the  various 
files  it  will  use  or  create  during  the  run. 


FILE 

f  or 003. dat 
f  or030. dat 
f  or015. dat 
f  or-090.  dat 


NAME 

e9vprb.scr 
e9vprb . dat 
e9vplt.dat 
escoplt. dat 


The  files  serve  the  following  purposes! 

egvprb. scr  —  contains  the  detailed  printed  output 
from  the  run.  Only  an  abbreviated 
version  is  sent  to  the  user's 
terminal  during  an  interactive  run. 
esvrrb.dat  —  input  data  for  egvprb,  containing  the 
f*  cylindrical  equilibrium  parameters. 

This  file  is  written  by  MHDEQUIL  for 
user-specified  equilibria.  A  separate 
routine  which  uses  the  output  of  a  rad- 
coupled  hydro  code  to  specify  the  equilibrium 
could  be  used  to  generate  this  file. 
f.  e9vplt.dat  —  input  data  for  the  plotting  code, 

EGVPLT ,  which  plots  the  answer  found  by 
egvprb.  This  file  is  generated  bv 
esvprb  when  it  is  used  in  the  mode  where 
the  code  finds  the  root, 
escorlt.dat —  input  data  for  the  plotting  code, 
fr  EGCOPLT ,  which  Plots  the  coefficients  aid 

trial  eigenvector  when  the  egvprb  code  is 
run  in  its  second  mode.  This  file  is  created 

■’  •'  ■"  •'  j.  A-' 


b y  esvprb  in  its  second  mode?  where  the  user 
manual  1y  searches  for  the  root. 

This  command  is  required  at  the  besinnins  of  each  session. 


STEP 


run  mhdequi 1 

This  code  sets  up  a  cYlindrical  equilibrium  for  esverb  based 
on  user  input.  All  data  is  entered  in  mks  units.  The  code  will 
prompt  the  user  for! 
awa  1 1 

btheta/radius/ifit 
pressure/ radius/ifit 
mass  densitY/radius/if it 

The  user  should  provide  free-format  inputs  consisting  of! 
line  1:  awall  —  wall  radius  (meters) 
line  2!  btheta  arraY  —  input  arraY  of  azimuthal 

masnetic  fields  (up  to  11  values) 
line  3:  radius  arraY  —  input  array  9ivin9  radii  at  which 

magnetic  fields  were  specified 
line  4!  ifit  array  —  input  array  of  0  or  1  for  specifying 

whether  linear  (ifit=0)  or  1/r 
(ifit*l>  interpolations 
are  to  be  used  to  specify 
the  magnetic  field  for  e9vprb. 
line  5!  pressure  ar-rav  —  pressure  profile  (nt/m**2) 
line  6:  radius  array  —  radii  where  pressure  specified 
line  7:  ifit  arraY  —  linear, 1/r  fit  switch  for  pressure 
line  8:  mass  density  —  density  profile  (kg/m#*3> 
line  9:  radius  arraY  —  radii  where  density  specified 
line  10:  ifit  array  —  linear, 1/r  switch  for  density 


line  4: 


radius  array 


line  5:  pressure  ar-rav  —  pressure  profile  (nt/m**2) 
line  6:  radius  array  —  radii  where  pressure  specified 
line  7:  ifit  array  —  linear, 1/r  fit  switch  for  pressure 
line  8:  mass  density  —  density  profile  (kg/m*#3> 
line  9:  radius  array  —  radii  where  density  specified 
®  line  10:  ifit  array  —  linear, 1/r  switch  for  density 

The  code  will  generate  input  data  for  esvprb,  and  store  it  in 
unformatted  form  on  esvprb.dat. 

This  step  is  needed  only  if  e9vprb.dat  does  not  already  exist. 

•  STEP. 3 

run  readSO 

This  code  reads  esvprb.dat  and  allows  the  user  to  see  the 
data  he  is  feeding  to  e3vprb. . 

9  STEP  4 

run  esvprb 

This  is  the  main  code.  It  will  read  the  equilibrium  data 
from  esvprb.dat,  and  will  prompt  the  user  for  additional 
information: 

r%  EVGUES, DEVAL, NEVAL, XK, XM, GAMMA, IBCL, IBCR,XL,XR. 

These  quantities  are  to  be  inputted  .in  free-format  as  a 
single  line  input.  They  stand  for  the  following  data: 

evsues  —  initial  9uess  for  the  eigenvalue,  normalized  as 
(omesatawal 1/vaa) #*2,  where  omega*#2  is  the 
eigenvalue  (a  squared  f requenc v ) , awa 1 1  is  the 
wall  radius,  and  vaa  is  the  Alfven  speed  at 
ths  wall. 

deval  —  eigenvalue  increment  for  use  when  the  code  is 
used  in  mode  2  —  the  code  will  examine  the 
system  for  neval  distinct  choices  of  the 
eigenvalue,  starting  with  evsues  and  incrementing 
£  the  eigenvalue  bv  deval  to  set  each  new  choice. 

neval  —  the  number  of  trial  eigenvalues  when  the  code 
is  used  in  mode  2.  If  neval=0  is  entered,  the 


code  will  run  in  the  first  mode,  disregarding 
d  e  v  a  1  and  using  e  v  9  u  e  s  as  t  h e  first  9 u e s s  for 
the  root  finder. 

xk  —  the  axial  wavenumber,  normalized  as  ktawall. 
xm  --  the  azimuthal  mode  number. 

samma  —  the  specific  heat  ratio,  default  value  is  b/3 
i  b  c  1  *  i  b  c  r  —  boundary  flags,  default  to  d 1  rich  let  bc's 
x  1 ,  x  r  —  9 rid  limits,  default  to  normalized  0,1  grid. 
The  code  will  also  prompt  (after  some  time)  for  a  parameter, 
SETNRM,  which  allows  the  user  to  alter  the  normalization  of  the 
determinant.  Typically,  this  parameter  will  be  specified  as  1, 
but  if  the  determinant  is  close  to  either  underflow  or  overflow 
on  the  VAX,  specifying  setnrm  different  from  1  may  allow  the 
calculation  to  proceed.  Alternatively,  the  number  of  nodes 
carried  by  code  can  t-e  reduced  to  avoid  overflow/underflow. 

The  output  from  the  code  will  be  on  egvprb.scr,  and  on 
plot  files  esvplt.dat  (mode  1)  or  egcoplt.dat  (mode  2). 

STEP  5 

run  e3coplt  (mode  2) 

No  user  input  is  required.  The  code  will  generate  plots  of 

1)  det  vs.  ev  —  the  charcteristic  determinant  vs 

the  eigenvalue. 

2)  a-coef f ic ient  vs.  r  for  each  eigenvalue  (3-d  plot) 

3)  b-c oef f i c i ent  vs.  r  for  each  eigenvalue  (3-d  plot) 

4)  c-c oef f i c i ent  vs.  r  for  each  eigenvalue  (3-d  plot) 

5)  eigenfunction  vs.  r  for  each  eigenvalue 

STEP  5' 

run  esvrlt  (mode  1) 

User  may  specify  which  plots  are  desired.  He  specifies: 
NGRAPHS  —  #  graphs  to  be  generated 

MGRAPH(J)  —  switch  for  each  type  of  graph 
1,0  mean  yes, no 

mgraph(l)  —  plot  a -coefficient  vs  r 
mgraph<2)  —  plot  b-c oef f i c i ent  vs  r 
mgraph(3)  —  plot  c-c oef f i c i ent  vs  r 
msraph(4)  —  plot  eigenfunction  vs  r 

* 

* 

* 

* 

% 


END  OF  EGVPRB. INF 


PROGRAM  MHDEQUIL 


C*** *#WRITES  MHD  EQUILIBRIUM  DATA 
C*****  FOR  EGVPRB 

C*****ON  FILE  F0R030.DAT 
C 

Farameter (np=21 7 

*  np2p-np+2f 

#  ninp=ll) 

d  intens  ion  bth  (nr2p  )  7  press  (  np2r  )  ■>  r  ho  (  hp2p  )  7 

#  va  (np2p)  »  cs  ( hp2p )  .  phvsrd  (iip2p)  » 

*  f  i np  < n inp )  ■»  r inp  ( ninp )  >  i f  i t  ( ni np ) 
data  awall/1. /, xmO/1 . 256637e-6/ 

data  gamma/ 1 . 66667/ 
data  x 1 >  xr/O. » 1 . / 
data  if i t/ninp*0/ 
r ew i n d  30 


1000 


print  1000 
f  ormat (lx» ’ enter 


tlx* 'enter  data  in  mks  units'/ 

3x . ' awal 1 ' / 

3x> ' b theta /radius/ if  it' / 

3x» 'pressure/radius/ifit'/ 

3xi ’mass  densi t Y/rad i us/ i f i t ' / 
lOx, ' ifit=(Oj 1 )  =  ( 1  inear  7 1/r)  fit 
>V)  awall 


data ' > 


read(5i#)  awall 
xr=awa 1 1 
del=xr-xl 
do  10  i  =  1 >  n p 

f ac= (float(i-l) /float (np-1 ) ) 
phvgrd ( i) =xl+del*f ac 
read(5i*>  finp 
read (57*)  rinp 
read(5j*>  ifit 

if  (rinp < 1 ) . ne.O. )  goto  500 

do  15  J  =2>  ninp 

if  (ri np  < J ) . eg. awal 1 )  rinp ( J )  =  1 . 0001 laws  1 1 

continue 

do  25  i  =  1 1 n p 

r=phY3rd(i) 
do  20  J=l»ninp 
if(rinp(J).le.r)  goto  20 

if<ifit<J)  .eg.O)  bth<i)**finp(j-l)  +  (finp(J)-finp(J-l)  ) 
*  *(r-rinp(j-l) )/ (rinp (J ) -rinp (j-1 ) ) 
if(ifit(J).eg.l)  bth  < i)  =f  ir»p  ( j-1 )  *rinp  < J-1)  /  r 
goto  25 
continue 
continue 


c***** 


do  30  i=l.ninp 
i  f i t ( i ) =0 
r i n  p ( i ) =0 . 
f i n  p ( i ) =0 . 
read (57*)  fine 
read  (5?  *)  r  ir  «• 
read(57»)  if: 
i f (r inp ( 1 ) . ne. 0 
do  32  J=l7ninp 
if(rinp(J).eg.awall) 
continue 
do  40  i*l7np 
r “phvsrd ( i ) 
do  35  J*l7ninp 


goto 


ri n  r ( J ) = 1 . 000 1  *  awa 1 1 


Ik 


* 


45 


47 


* 


50 


•  1001 

* 


•  1002 
60 


500 

»  1003 

* 

600 


iftrinp(J).le.r)  a  o  t  o  35 

i  f  (  i  f  i  t  ( j  )  .  e  a .  0 )  f  r  e  s  s  <  i  )  =  f  i  n  f  ( j  - 1 )  +  ( f  l  n  p  ( j  )  -  f  i  n  f  ( J  - 1 )  ) 

#  (  r-r  i  HF-  ( J  -1  >  )  /  ( r  i  n f  (j)-rinp(j-l)) 
if(ifit(J)  .  es  ,  1 )  press  (  i  )  =f  i rip  (J -1 )  #r  i  nr  ( j  -1 )  /  r 
3  o  t o  40 
continue 
continue 

do  45  i  =  l »  ninp 
i f  i  t  < i ) =0 
r  i  n  p  (  i  )  =0 . 
f i n  f  ( i ) =0 . 
read(5»*)  fine 
read (5»  *)  rinp 
read (5.*)  ifit 

i  f  ( r  i  n  p  ( 1 )  .  n  e .  0 .  )  soto  500 
do  47  J =2* ninp 

if ( r i ne ( J ) . ea . awa 1 1 )  r i nF  < J ) =1 . 0001  Hawaii 

continue 

do  55  i  =  l*riF- 

r- = f-  h  v  3  r-  d  ( i ) 

do  50  j  =  1 »  n i n p 

if(rinp<j).le.r)  3  o  t  o  50 

if(ifit(J).e«i.O)  rho  ( i )  inp  ( J-l )  +  (f  inp  ( j  )  -f  inp  ( J-l )  > 
#(r-rinp(J-l))/(rinp<J)-rinp<j-l)) 
if  (  if  it  (J  )  .  ea.  1 )  rho  <  i )  -f  inp  (J-l )  *r-inp  ( J-l )  /r 
3  o  t  o  55 
continue 
continue 
write (6> 1001 ) 

formatdhl* 19x, ' r ' , 12x, ’ bth' > 12x»  '  rho' » lOx, ' press'  » 

13x» '  va  '  > 13x. ' cs ' /) 
do  60  i  =  1 ?  n p 

va ( i ) =sar  t (bth ( i ) **2/xm0/rho ( i ) ) 
cs  (  i  )  =s«trt  ( samma# press  (  i  )  /rho  <  i  )  ) 

wri te  <6> 1002)  phY3rd (i) i bth(i) , rho (i) >Fress (i) i va (i) ics (i) 

format (5x  >  6el5. 5) 

continue 

write(30)  awa 11 »  bth »  rho»  pres s»  va»  cs 
endfile  30 
3  o  t  o  600 

write (6» 1003)  rinp(l) 

format  ( lhOi  5x»  '  #  #  *  #  #  #  #  #  #  *  input  errorl****r  inp  ( 1 )  =0.  is  expected' 
' ********** input  has  r inp ( 1 ) =' » el5. 5/) 
continue 


F-rosram  read  15 
parameter-  <n 3d p  =  101 ) 

dimension  p 1 1  s  r d ( n  s  d  p ) »pltr(nsdp) »  p 1 1 i  (nsdp) 
rewind  15 
do  100  m=i,4 

r  ead ( 15)  nsd»  pi t srd *  p 1 tr , f 1 1 i *  vmirn  vmax 
wri  te  <6i  1000)  ri3d 
f  or-mat  ( lhl  ?  5x *  ' n3d= ' > i6) 
wr i te (6> 1001 )  (pltsrd(i),  i=l»nsdp) 
format  <5x, ' p 1 tsrd= ' /50 ( lOx, 6el5.  5/ ) ) 
wr i te (6 » 1002)  <pltr(i)»  i=l»nsdp) 
f  or mat <5x> ' Pi tr  = ' /50 (10x,6el5.5/) ) 
wr  i  te  (6,  1003)  (plti(i),  i=l>nsdp) 
f  or  mat (5x, ' plti  =  ' /50 ( 10xi 6el5. 5/) ) 
wr i te (6> 1004)  Ymin, vmax 

f  orma t  <5x » ' ymi n= ' , el 5. 5> 5x > ' vmax= ' . el5. 5) 
continue 


p  r  o  9  r  a  m  read30 

%*%%%%%%%**.%*.***%****)■%*  * ##*##*#* 
c  reads  &  prints  for030  written  bv  * 

c  b y  mhde^uil  to  provide  data  to  esvrrb  # 

parameter < np=21 7 
#  np2p=np+2) 


10 

1000 

1001 


1002 

20 


dimension  bth (np2p) j  press (np2p) 7  rho (op2p) 7 

*  va (np2p) >  cs (np2p>  7  phY3rd (np2p> 
r  ew i n  d  30 

read (30)  awa 1 1 7  bt h 7  r ho  7  pr ess  7  va  7  c  s 
do  10  i=l7np 

fac= (float ( i-1 ) /float  (np-1 ) > 
phv3r  d  ( i ) =awa 1 1 *f ac 
wr i te (67 1000)  awa 11 

format ( lhO? 20x7 1  data  for  es vprb * / lOx 7 ' awa 1 1= ' , el5. 5/ ) 
wr i te (67 1001 ) 

format (5k 7 ' i ' 7  9x7 ' r ' 7  7x7 'bth* 7  7x7 ' rho ' 7 5x7 'press'  7 

*  Sx  7 ' va ' 7  8x  7 ' c  s ' ) 
do  20. i=l 7 np 

write (67 1002)  i 7  phvsrd ( i ) *  bth ( i ) 7rho(i) 7  press (i) 7va(i) 7  cs (i) 

f orma t ( 1 x  7 i5?  6el0. 3) 

continue 

•tint- 


PROGRAM  EGVPLT 
IMPLICIT  REAL #8 (D) 

REAL  1X1, 1X2, IY1, IY2 
I NTEGER  OTAPE ,  BUFFER ( 1 ) 

LOGICAL  LTIME 
PARAMETER ( NGDP= 101) 

D I MENS I ON  DPLTGRD < NGDP ) , DPLTR <  NGDP ) , DPLT I ( NGDP )  , 

*  PLTGRD (NGDP) , PLTR (NGDP) ,PLTI (NGDP) »  MGRAPH (4) 

INTEGER  T I TL I N ( 1 ) ,  T I TEND ( 1 ) 

LENBUF=1 
0TAPE=106 
NGRAPHS=1 
X 1R=150. 

X2R=750. 

Y1R=100. 

Y2R=700. 

N0VX=20 
NDVY=20 
LTIME=. FALSE. 

TIME=0. 

PRINT  900 

FORMAT (IX,' NGRAPHS , MGRAPH ' / 

*  5X , ' MGRAPH (J)=(1,0)=(Y,N) '  / 

*5X,'  J=1 ,2*3, 4= A, B, C, EV' ) 

READ (5,*)  NGRAPHS, MGRAPH 

CALL  GRAF I T ( 0 , OTAPE , BUFFER , T I TL I N ) 

CALL  GRAF I T ( 8 , OTAPE , BUFFER , - 1 ) 

REWIND  15 
NREAD=0 

DO  100  JGRAPH-1, NGRAPHS 
NREAD=NREAD+1 

READ  <15, END= 1 50 )  NGD , DPLTGRD , DPLTR , DPLT I , DYMIN, DYMAX 

I F ( MGRAPH ( NRE AD ) . EQ . 0 )  GOTO  5 

YMIN=DYMIN 

YMAX--DYMAX 

YYMAX=MAX (ABS(YMAX) , ABS(YMIN) ) 

IF ( YYMAX . EQ. 0. )  YYMAX=1 . 

YMAX=YMAX/YYMAX 
YMIN=YMIN/ YYMAX 
DO  10  J=1 , NGDP 
PLTGRD ( J) =DPLTGRD ( J) 

PLTR  <  J ) =DPLTR ( J ) 

PLTI <J)=DPLTI (J) 

DO  20  J=1 , NGDP 
PLTR  < J) =PLTR  < J) /YYMAX 
PLTI ( J) =PLTI (J) /YYMAX 
CONTINUE 
X 1=PLTDRD ( 1 ) 

X2=PLTGRD(NGD) 

IX1=X1 

I  X2=X2 

Yl-YMIN 

Y2»=YMAX 

IY1»Y1 

IY2«Y2 

CALL  GRIGEN ( X 1 , X2, X 1R, X2R, Y1 ,  Y2,  Y1R,  Y2R, 

*  '  R  '  LPMA  *. ' , 

*  TIME, I X 1 , 1X2, IY1, I Y2, NDVX , NDVY , LTIME) 

CALL  PLOTL (PLTGRD, PLTR, NGD, ’  $.  '  ) 

CALL  PLOTL (PLTGRD, PLTI, NGD, ’  ...*.') 

CALL  GRAF I T ( 4 , OTAPE , BUFFER , JGR APH ) 

CONTINUE 


PROGRAM  EGCOPLT 
c omp  lex  ev ,  det 
REAL  1X1, 1X2, IY1, IY2 
INTEGER  OTAPE, BUFFER ( 1 ) 

LOGICAL  LTIME 
PARAMETER  <  NGDP= 1 0 1 , 
t  nevalp=20) 

DIMENSION  wmin  (4) ,  YYmax  (4)  ,  detv  (neva  If)  ,  evv  ( neva  If)  , 

*  aa(4,nevalF,  n  3d  F-),xdata(  neval  f,  risdF)>Ydata(nevalp»nsdp) 

*  x  f  1  o  t  ( n  3  d  f  )  ,  y  f  1  o  t  ( n  3  d  F  )  , 

*  PLTGRD (NGDP) »  PLTR (NGDP) ,PLTI ( NGDP ) , MGRAPH (4) 

INTEGER  TITLIN(l) ,TI TEND ( 1 ) 

data  sph i / . 866/ , c f hi /-. 5/ 

LENBUF=1 
0TAPE=106 
NGRAPHS=1 
X 1R=150. 

X2R-750. 

Y1R=100. 

Y2R=700. 

NDVX=20 
NDVY=20 
LTIME=. FALSE. 

TIME=0. 

CALL  GRAF I T ( 0 , OTAPE , BUFFER , T I TL I N ) 

CALL  GRAF IT (8, OTAPE, BUFFER, -1) 

REWIND  90 
read (90)  neval 
d  o  8000  m=  1 , 4 
wmin  (m) =0. 

8000  YYmax (m) =0. 
detmin=0. 
detmax=0. 
do  8010  i=l, neval 
read (90)  ev»det 
evv ( i ) =real (ev) 
detv(i)=real (det) 
detmin=min (detmin,detv(i) ) 
detmax=max (detmax , detv ( i ) ) 
do  8008  M=1 , 4 

read (90)  nsd, Pi tsrd, p 1 tr , pit i , vmi n, Ymax 
do  8005  j=l,nsd 
8005  aa  (m,  i ,  J  )  =F-ltr  ( J  ) 

wmin  (m)  =min  (wmin  <m) ,  vmin) 

8008  YYmax (m) =max (YYmax (m) , Ymax) 

8010  continue 

do  8015  1*1=1 » 4 

Ymaxx=max  (at>s  (YYitiin  (m)  )  ,abs  (YYmax  (m)  )  ) 
wmin  (m)  =Yvmin  (m)  /  Ymaxx 
YYmax (m) = YYmax (m) /ymaxx 
do  8015  i=l, neval 
do  8015  J  =  1 , n s d 

8015  a  a ( m , i , J ) =a  a ( m , i , J ) / yma  x  x 
c 

ddmax=max (abs (detmin) , abs (detmax) ) 
do  20  i=l, neval 
20  detv ( i ) =detv ( i ) /ddmax 
detmax=detmax/ddmax 
detmi n=detmi n/ddmax 
evmin*min (evv(l) . evv (neval ) ) 
evmax=max ( evv ( 1 ) , evv ( neva 1 ) ) 
evmaxx=ma* (abs ( evmi n ) , at s (evmax ) > 


8020 


8050 

8100 


8200 


do  8020  i=l, neval 
evv ( i ) =evv ( x ) /evmaxx 
evmax=evmax/evmaxx 
e  vmi  n-evmi  n/evmaxx 
evrni  n=min  (  evmi  n  ,  0 .  ) 
evmax=rnax  (evmax.O.  ) 
detmin=min (detmin, 0. ) 
detmax=max  (detmax, 0.  ) 

X l=evmi n 
X2=evmax 
IX1=X1 
I X2=X2 
Yl=detmi n 
Y2=detmax 
IY1=Y1 
I Y2=Y2 
Jsraph=l 

CALL  GRIGEN(X1, X2, X1R,X2R, Yl, Y2, Y1R, Y2R 

*  ‘  EV  TED  *. ' , 

*  TIME, I XI. 1X2, I Yl, I Y2 , NDVX , NDVY , LT I ME ) 

CALL  PLOTL (evv, detv, neval , '  $. ' ) 

xf lot  ( 1 ) =evmin 
xf  1  ot (2) =evmax 
vpl ot ( 1 ) =0. 
yplot (2)=0. 

call  Plot  1  ( xF-lot ,  YFl  ot ,  2,  ' 

CALL  GR AF I T ( 4 , OT APE , BUFFER , JGRAPH ) 

CONTINUE 

do  8500  M=1 , 4 

Ymin=0. 

Ymax=0. 

xmin=0. 

xmax=0. 

do  8100  i=l, neval 
do  8050  J  =  l,ri9d 

xdata  < i , j ) =p It srd  U ) +evv ( i ) fcphi 
Ydata(i,J)=aa(m,i,J)-evv(i)*sphi 
xmin=min  (xmin,  xdata  (  i  ,  j  ) ) 
xmax=max (xmax, xdata ( i , J ) ) 
vmi n=mi n ( vmin, Ydata ( i , J ) ) 

Yinax=iiiax  ( vmax ,  Ydata  (  i  ,  J  )  ) 

continue 

continue 

xmaxx=max  (abs  <  xmi n )  ,  abs  ( xmax )  ) 

Ymaxx=max <abs< vmin), abs (vmax) ) 

do  8200  i=l, neval 

d  o  8200  J  =  1 , n  9  d 

xdata  < i , J ) -xdata ( i , J ) /xmaxx 

Ydata(i,J)=vdata(i>J) /vmaxx 

xmax -xmax /xmaxx 

xmi n=xm in /xmaxx 

Ymax=Ymax/ Ymaxx 

vmi n=Ymin/Ymaxx 

xmin-min (xmin,0. ) 

xmax=max <xmax,0. ) 

vm i n=m i n ( vmi n , 0. ) 

vmax=max ( vmax, 0. ) 

xl=xmin 

x2=xmax 

vl=Ymin 

v2*Ymax 

1 x l*x 1 


8250 


8300 


8400 


1 V 1=Y 1 
i y2-y2 

CALL  GRIGENtXl, X2, X1R, X2R, Y1 , Y2, Y1R, Y2R, 

*  '  R  LPMA  %.  '  , 

*  TIME, I XI, 1X2, I Yl, I Y2, NDVX  ?  NDVY ,  LT IME) 
do  8300  i=l?neval 

do  8250  j  =  1 , n 3 d 
xplot(j)=xdata(i?J) 
y  F- 1  o  t  <  J  )  =  y  d  a  t  a  ( i  ,  J  ) 

call  plotl(xplot?YPlot,n9d»'  $.') 
continue 

do  8400  i=l,neval 

xf-1 ot  < 1 ) =evv  < i ) #cf  h i /XMAXX 

XPl Ot (2) =XP 1 Ot (1) 

yp 1 ot ( 1 )  =  < YYmin (m) -evv ( i ) #sehi ) /YMAXX 

yf-  1  ot  (2)  =  (  y  Ymax  (m)  -ev v  (  i  >  #sph i  )  /YMAX  X 

x  f-  1  o  t  ( 1 )  =ma  x  ( x  p  1  o  t  ( 1 )  ,  xm  i  n ) 

xp  1  o t  (2)  =mi  n  <  xf-  1  ot  (2)  ,  xmax ) 

y  f-  1  o  t  <  1 )  =rna  x  ( v  f-  1  o  t  ( 1 )  7  yiti  i  n ) 

y  F- 1  o  t  ( 2 )  =m  i  n  ( v  p  1  o  t  ( 2 )  ,  y  ma  x ) 

call  F-lotl(xF-lot7YPlot,2,'  $ .  '  ) 

xplot  < 1 ) =ev v  (i)  3(c  c  F-  h  i  /XMAXX 

xplot<2)=<Fltsrd(n3d)+evv(i)#cphi) /XMAXX 

YF-lot  ( 1  )=-evv  (  i  )  *sphi /YMAXX 

YF- 1  ot  (2) =-evv  (  i  )  #sphi  /  YMAXX 

x  p 1 o  t ( 1 ) =ma  x ( x  p 1 o  t ( 1 ) , xm i n ) 

xplot  (2)=min (xplot (2) 7  xmax) 

y  p 1 0 1 ( 1 ) =ma  x(YPlot<l)»Ymin) 

YPlot(2)=min(YPlot<2)7Ymax) 

call  plotl  (xplot,  YF-lot,  2.  '  *.  '  ) 

continue 

xf  1  ot  ( 1 ) =evmin*c phi /XMAXX 

xplot (2) =evmax*cPhi/XMAXX 

vp lot ( 1 )=-evmin#sphi/ YMAXX 

yplot  (2) =-evmax*SPhi /YMAXX 

call  f  loti  (xplot,  YF-lot,  2,  '  *.  '  > 

CALL  PLOTC (0. , 0. , 1 , ' 

XPLOT ( 1 ) =EVMIN*CPHI 

YPLOT ( 1 ) =YYMAX (M) -EVMIN*SPHI 

X PLOT ( 2 ) =E VMA X # CPH I 

YPLOT (2) =YYMAX (M) -EVMAX*SPHI 

XPLOT  (3)  =PLTGRD  (NGD)  +EVMAX*CPHI 

YPLOT (3)=YPLQT (2) 

XPLOT (4) =PLTGRD (NGD) +EVMIN#CPHI 
YPLOT ( 4 ) =YPLOT ( 1 ) 

XPLOT (5)=XPL0T(1) 

YPLOT (5) = YPLOT ( 1 ) 

XPLOT  (6)=XF’L0T  (5) 

YPLOT  (6)  =YYMIN (M)  -EVMIN*SPHI 
XPLOT (7)=XPL0T (4) 

YPLOT ( 7 ) =YPLOT ( 6 ) 

XPLOT  (8)  =XF'LOT  (3) 

YPLOT ( 8  >  =YYM I N ( M ) -E VMAX  #SPH I 
XPLOT (9)=XPL0T (2) 

YPLOT ( 9 ) =YPLOT  <  8  > 

XPLOT  (10)  *=XPLOT  (6) 

YPLOT (10) = YPLOT  <  6 ) 

XPLOT (11) =XPLOT (7) 

YPLOT (11) =YPLOT ( 7 ) 

'  XPLOT ( 12) «*XPLOT (4) 

YPLOT (12) “YPLOT  <  4 ) 

XPLOT (13) =XPLOT (3) 

YPLOT (13) =YPLOT (3) 


ft 


\  vO  \  v  v'  •;  '/* .  '  -  «•  ».* 


XF'LOT  (14)  =XF'LOT  (8) 

YPLOT  (14)  =YF'LOT  (S) 

XF'LOT  ( 15)  =XF‘LOT  (9) 

YF'LOT  (15)  =YF’LOT  ( 9 ) 

XPLOT (16) =XPLOT (6) 

YPLOT  (16)  = YPLOT  (6) 

DO  8900  11=1, 16 

IF  (XPLOT  (II)  .GT.  XMAX)  XF'LOT  ( 1 1  >  =XMAX 
IF ( XPLOT (II).LT.XMIN)  XPLOT ( 1 1 ) =XMIN 
IF  (YF'LOT  (II)  .  GT.  YMAX )  YPLOT < 1 1 ) =YMAX 
IF (YPLOT (II) . LT. YMIN)  YPLOT < 1 1 ) =YMIN 
CONTINUE 

CALL  PLOTL (XPLOT, YPLOT, 16, *  ...$.') 

J  sraph=ni+l 

call  9rafit(4,otape,buffer,JsraPh) 
continue 

call  srafit(9,o tape, buffer, tit end) 

stop 

end 


c ########################################################### 
c 

subroutine  bc(ibndrv) 

C  IMPLICIT  REAL*8(A-H,Q-Z) 


PARAMETER (NP  =  21, 

3  nesF  =  1, 

3  nsir  =  1, 

)  nd i p=  1 , 

5  ntip=  4, 

j  03dF=101, 

7  nm2p=np-2,  nml p  =of-1  ,  opI  f-=of  +  1 ,  of  2f  =op+2, 

7  nesssF^nesF-tnesp,  neF  =4*ne<?p+l, 

3  maxp=nesp*03dp,  nsap=16*np2p,  oa  1  p=7*oe*»5‘»p 

)  odrti2p=ne3p*OF-2p,  ndm3p=ne3S3p*np2p , 

5  odfii23p=3*oe«tP ,  ndm33p=3*ne3SSF) 


common/baksub/aa  <ndm2p, odm2p) 

commoo/esnvcs/esvs (op2p, nesp) , usi (ndm2p)  ,  wa (ndm2p) 
c ommon/bdYc ds/bc 1 (3, nepp) , bcr  <3, oppp) 


common/matrix/ 


corrirrion/tosf  rm/ 


a ( of  2p , nepp, nepp) , 
b  (op2p,  oe«pp,  ne«»p) , 
c (np2p,nepp,nepp) 
sp ( op2p , oep p , oepp ) , 
br  (op2p,  nesp,  nesp) , 
cr  (of-2p,  omp,  nepp) , 


ai  (np2p,  oepp,  nesp)  , 
bi <op2p, nepp,nepp) , 
c i (op2p, nepp.nepp), 


3  ars(np2p,nepp,riepp)  ,ais  (of-2p»  nepp.nepp), 

l  brs  (op2p,  oepp,  ne«pp) ,  bis  (np2p,  nepp,  nepp) , 

>  crs  (op2p,  ne«tp,  nepp)  ,  c  is  (op2p,  oepp,  nepp) 

commoo/5F  vals/F  hi  (maxp)  ,  spvals  (nadp,  of-2p) 
c  ommo  o  /  9rids/n9d,xl«xr,phv9rd(np>*pltsrd(n9dp) 
commoo/intser/nep, oppsp, ndm2, ndmS, ndm23, ndm33,max, itrmax 

COMPLEX  ev,evold,esvs,usi,cmxe,aa,al,phi 
COMPLEX  a,  b,  c  ,  det,  detnrm 
COMPLEX  be  1 , bcr, z norm 


common/nc  l/nrri2,  nml ,  o,  opI,  of  2 
c  ommo  n/rcl/r(np),rn(op) 
common/sac l/sa (nsap) 

commoo/phiOl /pOI (op) /rhi02/p02(np) /phi03/p03(np) 
common /phi 1 l/pll(np)/phil2/pl2(op)/phil3/pl3(op) 
commoo/phi21/p21 (of) /phi22/p22(op) /phi23/p23(np) 
c  ommon/phi i l/pil(op)/phii2/pi2(np) 
c  ommo o /phi i3/pi3(np)/phii4/pi4(op) 
c  ommoo/ph i4/e ( op2p ) /phi5/f (op2p ) 
commoo/bodv  1  s/bcOp,  be  Of-1 ,  be  1  p  1 ,  be  Ip,  be  1 ,  bcO 
commoo/errc 1/err (op2f) 
commoo/serrc 1 /amxer 1 , amxer2, ermax 
common /duml /dl ( op2p ) /dum2/d2 (  op2f-  ) 


c 

c#### 

c#### 

c#### 

c#### 


commons  for  spline  iote9rals 

ssi (op2, i) , dsi ( op 2, 7, J ) , tsi (of2,49, k ) 
i^oumber  of  siosle  iotesrals 


J  = number 
k=number 


double 

triple 


iotesrals 
i ote9r a  1  s 


common/ssic  1/ssi  (np2p,  nsip) 
c  ommo n  /dsicl/dsi  ( n  f-2f  ,  7  ,  n d  i  p  ) 
common/tsi c 1 /ts i (np2p, 49, nt ip) 

c  ommon/nss i f /ns i , nsJ (nsiF) , ns v (2, ns i f ) 
common/ndsip/ndi,ndJ  ( nd  i  p )  , ndv (4, nd iF ) 
c  orrirtion/nts  i  p/nt i , nt J  (nt iF )  ,  ntv  (6,  nt  i  f  ) 

equivalence  (sa,sa3) 
dimension  sa3 (4, 4, ne2e) 


setup  boundary  conditions 

ibndr-Y  =  1  left  boundary  condition 
i bn dry  =  2  risht  boundary  condition 

COMPLEX  tembc (3) , srf trm 
index  =  ( ibndrv-1 ) Inpllnei 

del i  =  1 . / (xr-xl ) 

do  90  i=l , neq 

90  to  (10,20)  ibndry 

10  continue 

do  11  j - 1 , 3 

tembc  < J )  =  be  1  ( j , i ) 

11  continue 
90  to  30 

20  continue 

do  21  J  =  1 , 3 

tembc ( J )  =  ber  ( J , i > 

21  continue 
30  continue 

tl  =  ABS (tembc (1) ) 
t2  =  ABS (tembc (2) ) 
t3  =  ABS (tembc (3) ) 

if (tl . ne. 0. . or. t2. ne. 0. )  so  to  40 
print  100 
WRITE (3, 100) 
stop  300 

40  if (tl. ne.O. . or. t3. ne.O. )  so  to  41 
so  to  90 

41  j f (t2. ne. 0. . or . t3. ne. 0. )  so  to  42 
do  50  k=l»ndm2 

aa  (  i ndex+i , k )  =  (0.,0.) 

50  continue 

aa ( i ndex+i , index+i )  =  (l.,0.) 
so  to  90 

42  if (tl. ne.O.)  so  to  43 
print  101 

WRITE (3, 101 ) 
stop  301 


43  if (t2. ne.O.)  so  to  44 


print  102 
WRITE (3, 102) 
stop  302 


44  if(t3.ne.O.)  so  to  45 
so  to  (61,62)  ibndry 
6 1  continue 


srftrm  =  -deli*a(  1  .,  i  ,  i  )  #bc0##3*tembc  ( 1 ) /tembc  (2) 
so  to  63 
62  continue 


srftrm  =  del  i  *'a  ( np2,  i  ,  i  )  #bc  1  *#3#tembc  ( 1 ) /tembc  (2) 


c 


c 

c 


63  continue 

aa ( i ndex+i , i ndex+i )  =  aa ( index+i , index+i )  +  srftrm 
so  to  90 

45  continue 
print  104 
WRITE (3, 104) 
stop  304 


90 

continue 

100 

format ( lx, ' 

stop 

300 

improper- 

boundary  conditions 

imposed '  ) 

101 

format (lx,' 

stop 

301 

bound  a  r  y 

conditions 

not 

yet 

implemented '  ) 

102 

format  < lx, ' 

stop 

302 

boundary 

conditions 

not 

yet 

imp  1 emented  '  ) 

103 

format ( lx, ' 

stop 

303 

boundary 

conditions 

not 

yet 

implemented '  ) 

104 

format  < 1 x , ' 

stop 

304 

boundary 

conditions 

not 

yet 

imp  1 emented ' ) 

c 

return 

end 

COMPLEX  function  c evalf (evtr i 1 , f cn) 

C  IMPLICIT  REAL*8(A-H,0-Z) 

c###*#  written  bv  J.  c.  macmahon  univ.  of  texas  at  austin  oct  1973 

c###]K#modif  ied  by  w.  h.  miner  univ.  of  texas  at  austin  Jun  1976 

c##**#modif  ied  by  a.  a.  mondelli  science  bppI  ications,  inc.  mar  1981 

c*##****##**  Oct  26  —  2.9 
c#¥*#**#*#*  oct  23  —  mate  with  fen  rks 


c 

external  fen 
c 

c  ommon/c evl f 1 /dxl ,ftest,dftest,dxtstl,dxtst2,tstep,tsuad 
common/c  evlf 2/ imax  , ieval ,  key , i s t ep  «  j start 
c  ommon/c  ev lf3/fl,xl 
c 

COMPLEX  dx  1 ,  x«t  (2)  ,  evtr  i  1 ,  f  cn 
COMPLEX  xO, x 1 , x2 
COMPLEX  f0,fl,f2 
COMPLEX  a0,al,a2 
COMPLEX  a,b,c,d,df,x 
c 

data  dxl7ftest7dftest7dxtstl7dxtst27tsteP7tsuad7 imax/ 

1  (.001,0.).  1.  c—08,  1 .  e-09,  1.  e-08,  1 .  ,  .  2,  .  02,  15/ 

c 
c 

c ################################################################### 
c 

c  dxl  initial  step 

c  ftest  tolerance  on  fen 

c  dftest  minimum  df  (Jexit=5) 

c  dxtstl  allowed  estimated  error  in  ev 

c  dxtst2  maximum  dx  (Jexit=4> 

c  tstep  defines  'small'  step 

c  tsuad  test  for  neishtorin?  root 


\  %  \  *.  ■  *  • 


imax  maximum  number  of  calls  to  fen 


c ################################################################### 


root  solver  modified  to  restrict  search  to  REAL  ROOTS  ONLY 


data  epsi2/  l.e-16  / 
data  Jstart/O/ 


DX1R=MIN (REAL (DX1 ) , ABS (REAL (EVTRIL) / 10. ) ) 

TQUAD=DX1R**2 

DX 1=CMPLX (DX 1R»  0. ) 

DXTST 1 =DX 1 R*  *2/ 1 000000. 
DXTST2=DX1R#*2*100. 

DFTEST=FTEST/DX1R**2 


WRITE  <3* 1 ) 

1  format (/'  i,  evr,evi,  abs9<fcn)»  kode'/) 


*********56^  uf  for  first  iteration 
2  istep=-l 
x2=evtri 1 
i eva 1=0 
key=0 
s  o  t  o  1 0 


***#***test  i eva 1 i  dx 
5  i f ( i eva 1 . 9e. imax ) 9o  to  99 
if (istep. e9.0) 90  to  6 


d=x-x2 

t=REAL(d) *#2  +  AIMAG<d)**2 
if  (t. It. dxtstl )  go  to  110 
if <t. 3t. dxtst2)  90  to  120 
*******shift  previous  values 
6  fO=f 1 
f  l=f2 
xO=xl 
x  1 =x2 

****new  x-values 
x2=x 

*****new  value  of  fen 
10  ieval=ieval+l 
f2=fcr.<x2) 

********  test  for  convergence 
af 2s=REAL (f 2) **2  +  AIMAG ( f 2) * *2 


WRITE (3, 11) ieval, x2,af2s,  istep 
1 l  format (Ixji3*2el3.5.5x7el0.2>lx»i4) 


14  if <af2s. It. f test)  go  to  100 

*****  test  for  mode  of  next  step 
if <istee)20,40» 15 


15  if  (t. It. tstep)  so  to  40 

*********  first  stePi  or  previous  step  large.  take  new  step  dxl 
20  if  (Jstart. e9. 1 ) 90  to  25 
x*x2+dx 1 
istep=0 
so  to  5 


.*  V  *  •  ■  •  *  •  V*  -  ’  *  '  *  *  *  *  ■  _“  •  -  *  9  *  -  *  ■  •  *  *  *  ‘  >  *  *  “  *  - 


’  -  •>*.  V- 


vJ 


J start=0 
i step=0 


40  df =  <  f 2-f 1 ) / ( x2-x 1 ) 

41  t=REAL<df ) **2  +AIMAG <  df ) *#2 
if (t. It. dftest)  30  to  130 

45  x=.5*<xl+x2-(fl+f2)/df) 
if(ist eF.se. Dso  to  50 
istep=istep+l 
so  to  5 

Hi  *  t  *  *  #  t  quadratic  extrapolation  from  points  x0ixl*x2 

50  aO=fO/ ( (xO-xl ) * <x0-x2)  ) 
al=f  1/  (  (xl-x2)  #  ( x  1-xO)  ) 
a2=f 2/ < (x2-x0) * (x2-xl ) ) 


a=aO+a 1 +a2 
i step  =0 

t=REAL(a)#*2  +  AIMAG<a)#*2 
if <t. It. epsi2) so  to  40 

b=aO* (xl+x2)  +al*(x2+x0)  +a2*(x0+xl> 
b=-b 

c=a0*xl*x2  +al*x2*x0  +a2#x0*xl 

d=SGRT ( b * *2-4 . O# a * c ) 

xq<l)=<-b+d)/<2.*a) 
xs (2) = (-b-d) / (2. *a ) 
xirl=real ( xq ( 1 ) ) 
xqr-2=r ea  1  ( xq  (2)  ) 
xq(l)=cffiplx(xqrl,0.  ) 
xq(2)=scfiiplx(xqr2>0.  ) 

d=xs ( 1 ) -x 

t=REAL(d)*#2  +  AIMAG(d>#*2 
d=xq (2) -x 

tl=REAL(d) *#2  +  A I MAG (d) ##2 
i  =  l 

if  <tl. st. t) so  to  55 

i=2 

t=tl 

55  df=2. *a*xq <i)  +b 
x=x  q ( i ) 
j.  step=l  +  i 
d=x q < 1 ) -xq  <2) 

tl=REAL(d)##2  +  AIMAG<d>**2 
if  (tl . st. tsuad) so  to  60 

WRITE <6i 66)  xq (3-i ) 

WRITE (3»  66)  XQ(3-I) 
i step=3+i 

60  if (t. It. dxtstl )  so  to  110 
30  to  5 


66  format ( 1 x» ’ war ni 09 »  neishborina  root  at  ' »  2  e  1 3 . 5  /  > 


c 

c 


exit 


O 


c 

c 


c 


m 


c 

#  C 

c 

c 


e 


******  ieval  .  st.  imax 
99  kev=4 
J ex i t=6 

* * * * * *  c  o river  sene  e 
100  cevalf=x2 

wr 1 te (6,  112)  ieval-,  x2,  i  step ,  ke  y  ,  J  ex  i  t 
wr  i t  e  <3, 112)  i eva 1 , x2, i step  ,  k  ev , J ex i t 
112  f or mat ( 1 x , i3, '  conversed  root= ’ , 2el3. 5* 

.  ’  i step , key , J exi t=  ',3i5/> 

r  e  t  u  r  n 

******  problems 

110  kev=l 

if  (  i  step.  1 1.  1  >.so  to  200 
c  eval f=x 

WRITE (6, 1 11) ieval,x, istee 
WRITE (3, 111)  IEVAL, X, ISTEP 
return 

111  format (lx, i3,2el3.5» '  del  x  .It.  dxtstl',i4) 

120  kev=2 

so  to  200 
130  k ev=3 

so  to  200 
140  kev=5 

so  to  200 

200  WRITE (6, 201) kev 
WRITE (3, 201)  KEY 
Jexit=kev+2 

if (key. es. 5)  so  to  5 
so  to  100 

201  format (lx,'  problem  kev=',i3) 
end 

COMPLEX  function  fcn(ev) 

IMPLICIT  REAL*8 ( A-H, 0-Z ) 


PARAMETER (NP  =21, 

2  nesp=  1, 

3  nsiF=  1, 

4  ndip=  1, 

5  ntip*  4, 

6  n3dp=101, 

7  nm2p=np-2, nml p=ne-l , np lF=np+l , np2p=np+2, 

7  nesssp=nesp*nese, nep=4*nesp+l , 

3  maxp=nesp*nsdp>nsap=16*np2p,nal p=7*nesssp, 

4  ndm2p=nesp*np2p»  ndm3p =nesssp*np2p , 

5  ndm23p=3*nesp, ndm33p=3*nessap) 


i 

i 


i 


* 


common/baksub/aa (ndm2p, ndm2p) 

common/esnvc  s/esvs (ne2p, nesp) , usi <ndm2F ) , wa (ndm2e) 
common/bdycds/bc I (3, nesp) ,bcr(3,nesp) 
c  ommon/ma tr i  x  /  a(np2p,nesp,nesp), 

1  b <np2p, nesp, nesp) , 

2  c (np2e, nesp , nesp) 

common/tnsf rm/  ar (nr2p, nesp, neap) ,  ai (np2e,  neap,  neap)  , 

1  br (np2p» neap, neap) ,  bi (np2p, neap, neap) , 

2  cr (np2p» neap , neap )  ,  c i (np2p, neap, near ) , 

3  ars (np2e , neap, neae ) , a  is  <np2p , neap , neap) , 

4  brs  (riF-2p,  rear,  neap)  ,  bis  (np2p,  neap,  neap)  , 


i  c r  s ( nF 2p 7  neap 7  neqp) 7  c 1 s ( np2p 7  nesp . neqp ) 

common/spvals/phi  (maxp)  7  srvals  (nsdFi  hf2f  ) 
common/srids/risd.  xl>  xr>  F-hvard  (np)  i  pltsrd  (nsdp) 
c  ommon/  i  nt  ser/nes ,  neFsi 7  ndm27  rid  m3  7  ndm23i  ndm33i  max  >  i  trmax 
COMMON/ESTORE/EVSTORE 

COMPLEX  ev>evoldie9vs,usi.cmxeiaai3liphi 
COMPLEX  a.b.c.detidetnriTi 
COMPLEX  bcl.bcr.znomi, EVSTORE 


c  ommo n /  n c  1  / nm2 >  nm  1  >  n  7  n p  1 7  n p2 
c  ommon/ re  l/r(np)  7  r-n  (np) 
common/sacl/sa  (nsap) 

c  ommori/F  hiOl  /  f-01  (np)  /phi02/p02  (np)  /p  hj.03/p03  ( np) 
c  ommon/F  hi  1 1  /p  1 1  ( np )  /ph i  12/p  12  ( np  >  /F  h  i  13/p  13  ( np  ) 
c  ommon/phi21/p21 (np) /phi22/p22 (np) /phi23/p23 (np) 
c ommon/ f- hi  i  1/pi  1  (np)  /  phi  i2/pi2  (np) 
c  ommon/phi i3/pi3 (np) /phi i4/pi4 (np) 
common/phi4/e(np2p) /phi5/f (np2p) 
c ommon/bndvls/hcOP7  bc0pl*bclpl»bclp>bcl*bc0 
common/errc  1/err  (np2p) 
c  ommon/serrc  1  /amxer  1 7  amxer-27  ermax 
common/duml/dl (np2p) /dum2/d2  <np2p) 
c 

c  commons  for-  spline  integrals 

c 

c####  ss  i  (riF  27  i  )  7  dsi  ( np27  7t  J  )  7  tsi  ( ne27  49?  1: ) 
c####  i=number  of  sin9le  integrals 
c####  J=number  of  double  integrals 
c####  k=number  of  triple  integrals 
c 

common/ssic 1 /ssi (np2p  7  nsip) 
common/dsicl/dsi(np2P777ndip) 
c  ommon /tsicl/tsi (np2P7  497  ntip) 
c 

common/nssip/nsi7  ns  j (nsip) 7  ns v (2?  nsip) 
common/nd si  p/nd i 7  ndJ (ndip) 7  ndv (4t  ndip) 
common/ntsip/rit  i  7  nt  J  (ntip)  7  ntv  (67  ntip) 
c 

equivalence  <sa7sa3) 
dimension  sa3 (4? 47 np2p ) 
c 

common  / los i  c  3. /I i ntal  7  ldtnriri7  lrerdr- 
logical  1 i ntal 7 ldtnrm? lrerdr 
c 

dimension  f-  r  n  1 1  ( 4 ) 

data  F-rntl/'  coe'7'ff  m'  7  '  atri  '  7  '  x  '/ 

COMPLEX  temp 
c 

c ommon /time/t 
data  immax/1/ 
c 

i f  ( . not. 1 intal )  90  to  31 
c 

c - initialize  values. 

c 

T=GETIME ( dum) 
c 

call  matrix (ev) 
c 

T=GETIME (dum) -t 


print  904  .t 
WRITE (3, 904)  T 


del-xr-x 1 
del i=l . /del 
del2i=deli*deli 
amxerr*.  05 
amxer  2=. 4 


r r  int  1 00 , x 1 7  x  r 
WRITE (3> 100)  XL »  XR 

100  f  ormat  ( lx,  '  values  initialized  xl=*tel2.5» 

— chose  9 rid. 

call  9 r i d  ( r » 1 ) 

-  setup  grid  for  plots  - 

do  50  i  =  1 »  n 9 d 

Plt9rd(i)=(i-1.0)/(n9d-1.0> 

50  continue 
print  102 
WRITE (3» 102) 

102  format (lx» ' 9rid  chosen’) 
print  101 , ( r ( i ) , i=l , n ) 

WRITE  (3»  101 )  ( R  <  I )  »  I  =»  1  *  N ) 

101  format (lx> el5. 7) 

— set  up  splines, 
call  bsplcf 

do  950  i«l,4 
do  950  J  =  1 1 4 

WR I TE ( 3 i 9 1 0 )  (sa3(i,J»k) ,H=l,ne2) 

950  continue 

T=GETIME ( dum) -t 
print  901  »t 
WRITE (3, 901)  T 

— compute  integrals, 
call  opset 

T=GETIME (dum) -t 
print  902»t 
WRITE (3, 902)  T 

— set  up  decomposer, 
call  deset(O) 

WRITE (3» 198) 

WRITE (3, 199)  bcOtbcOp, bcOpl.bclPl, bclp.bcl 
198  format (lx«'  spline  boundary  values') 

T=GETIME ( dum) -t 
print  903. t 
WRITE (3»  903)  T 

31  continue 


temp*  < 1 . 1 0. ) 


initialise  ei sen vectors 


do  92  1  =  1 1  n d Hi 2 
usi(i)  = ( 1 . i 0. ) 

92  continue 

r  r  i  n  t  1 05 
WRITE <3, 105) 

105  f ormat ( 1 x» *  ini tia 1  eigenvectors' ) 
WRITE (3» 104)  <  US I ( I ) » I  =  1 » NDM2 ) 

94  continue 

104  format<lx,2el2.5) 


- initialize  error  vector - 

do  93  i  =  l,n 
err ( i ) =0. 

93  continue 

WR I TE (  3  ■>  106) 

106  f  ormat  C 1  x  >'  error-  vector') 
WRITE  <  3i  101 )  ( err  <  i )  >  i=*l »  n> 

-  calculate  matrix  elements 

im=l 

83  continue 


-  store  physical  mesh 

do  42  i  =  1  •>  n 
rhv9rd(i)=r(i)#del+xl 
42  continue 


call  matrix (ev) 

T=GETIME ( dum) -t 
f-  r  i  n  t  904 ,  t 
WRITE (3»  904)  T 

do  11  i - 1 »  n 
do  11  J  =  1  >  n  e  s 
do  11  k=lineg 
ar  < i . j , k ) =  REAL (a ( i ,  J  ,  k  >  > 
ai <i,J,k)=AIMAG<a<i, J,k> ) 
br  <i, J  ,  k>*  REAL (b <i, J , k) ) 
t-i  (i,  J  ,  k)=AIMAG<b  <i,  J  ,  k)  ) 
c  r  < i , J , k ) -  REAL ( c ( i >  J . k  > ) 
c  i  <  i  ,  J  »  k  )  “AIliAG  ( c  ( i ,  J  ,  k )  ) 

11  continue 


-  setur  matrix  elements  in  selines 

do  22  j=l i nes 
do  21  k=l . nei 

call  depse(ar(l»iik)iars(l*iik)) 
call  derse(ai(l.i)k).ais(l.iik)) 
call  derse  ( t>r  1. 1 ,  i  i  k  )  i  brs  ( 1 » i ;  k ) ) 


vv- .  •  .*• 

s,  .  .  ‘  .  V  •**  • 


call  depse(bi(l»i»k)»bis(l»i.k)) 
call  depse(cr(liiik)icrs(l,iik)) 
call  deFse(ci(l.iik)»cis(l,i,k)) 

21  continue 

22  continue 

do  190  k  =  l.ne«i 
do  190  J  =  1 7  n  e  s 
do  191  i  =  1 7  n 

WRITE (3, 199)  ar  <i,J,k> »ai  (i,J,k) 
1  b  r  <  i  .  J  .  k  )  ,  b  i  <  i  7  J  .  k  )  > 

1  c  r  (  i  7  J  7  k  )  7  c  i  ( i  7  J  7  k  ) 

191  continue 

do  192  i  =  l  7  riF  2 

WRITE  (3,  199)  ars  <i,  J  ,  k)  .ais  (i,  J  , 
1  t«  r  s  <  i  7  J  >  k  )  i  b  i  s  (  i  >  J  >  k  )  7 

1  c  -  s  <  i  7  J  7  k )  7  c  i  s  ( i  7  J  7  k ) 

192  continue 
190  continue 

199  f  or-mat  ( 1  x  7  6el2.  5) 


T =GET I ME ( d  um ) - 1 
f  r  i  nt  905?  t 
WRITE (3, 905)  T 

if < . not. 1 intal )  so  to  81 

i-f  (.  not.  lrerdr)  90  to  85 

- calculate  new  9rid  if  necessary 

if < im. se. immax)  so  to  85 

do  82  i=l7nes 

call  sel err (ars < 1 7 i 7 i ) 7  err) 
call  sf lerr (a i s ( 1 7 i 7 i ) 7  err ) 
call  selerr (brs (1 7 i 7 i ) 7  err ) 
call  sf 1  err ( bi s ( 1 7 i 7 i ) 7  err ) 
call  sf-1  er r  (crs  ( 1 7  i  7  i )  7  err ) 
call  selerr (c is ( 1 7 i 7 i ) 7  err ) 

82  continue 

WRITE <3i 6000)  (err(mm)  ,mm=l,n) 
i sr id=2 

if(im.ne.O)  isrid=3 
call  srid  (err-7  isrid) 
call  rmove(err) 
call  bsrlcf 
call  deset(O) 
im=im+l 

WRITE  (3>  6000)  <  r  (mm)  ,  mm=l ,  n ) 
6000  f ormat ( 1 x 7 10el2. 3) 

if (im. It. immax)  so  to  83 

- end  of  knot  adjustment - 

85  continue 

- calculate  inte9rals - 


call  sfeval 


T=GETIME (dum) ~t 
print  906 > t 
WRITE (3i 906)  T 


c 

c 

c  setup  spline  values  for-  reconstruction 

do  41  j=l*np2 

call  sPvl(j7Pltsrd<l)iSPvals(l»J)iri9d) 
41  continue 
c 

lintal=. false, 
c 

81  continue 

c 

c - 

c 

do  27  l;=l,ne«i 
do  27  J  =  l,ne«i 
do  27  i  =  1 , n p 2 

a ( i 7  J  7  k ) =DCMPLX (ar s ( i 7  J  7  k ) 7  a i s  ( i  7  J 7  k ) > 
t.  (  i  ,  J  ,  k  )  =DCMPLX  ( br s  < i , J  7  k  > , b i s  < i , J , k  >  > 
c ( i 7  J  7  k ) =DCMPLX ( c  r  s  (  i  7  J  7  k  )  7  c  i  s  <  i  7  J  7  k  >  > 
27  continue 
c 

call  setbc  <ev) 


c  zero  coefficient  matrix  before  each  iteration 

c 

do  70  i=l7ndm2 
do  70  .i  =  l7ndm2 
aa ( i 7  J  )  =  (0. 7  0.  ) 

70  continue 
c 
c 

do  23  k  =  1 , n p 2 
c 

imin=max0(k-37 1) 
imax=min0<k+37  nr2> 
do  26  i=imin7imax 
i l=i-k+4 

.imin-maxO  ( 1 7  k-3) 

Jmax=minO (np27  k+3) 
do  20  1  =  1 7ne«» 
do  20  m= 1 7  n e 9 
do  24  j=Jmiri7  Jmax 
J 1 = j -k  +4 
i2=il+7*  1-1) 


i  n  d  r  =  ( k  - 1 )  #  n  e s  +  1 
i n  d  c  =  (i~l)#nes  +  m 
aa(indr> indc)  =  aa(indr.indc) 

1  -  d  e  1 2  i  #  a  ( J  ,  1 ,  m )  *  ( t  s  i  ( k  ,  1 2 ,  1  >  + 1  s  i  (  k  ,  i 

2  +  deli*b (J 7 l,m) *  tsi(k,i2,3) 

3  +  c(J.  l,m>*  tsi  (k, 12.4) 

24  continue 
20  continue 
26  continue 

setup  boundary  conditions 

i  f  (  k  .  n e .  1 .  a n d .  k  .  n e .  n p 2 )  so  to  71 
if (k . ne. 1 )  so  to  72 
call  t<c(l) 

30  to  71 
72  continue 
call  be (2) 

71  continue 
23  continue 

if (. not. ldtnrm)  call  setnrm 
if (ldtnrm)  so  to  970 

call  F-rnt  (pr-nt  1 ,  aa,  2#ndm2,  ndm2,  1 . 2#ndm2,  1 ,  ndm2) 

970  continue 

call  matnrm 

call  le*itlc(aa-7  ndm2,  ndm27  esvs,  1 7  ndm27  1 .  wa .  det .  ier) 
if (ldtnrm)  90  to  980 

call  ft  nt  <f  rntl 7  aa> 2#ndm2, ndm2. 1 . 2*ndm2, 1 7  ndm2) 

980  continue 

temp  =  det 


normalization  of  determinant 

if (. not . ldtnrm)  90  to  90 
t  emp=t  emp / d  et  nrm 
90  to  91 

90  continue 
detnrm=temp 
ldtnrm*. true. 

91  continue 
nevas=2 

print  700  7  e v . t  emp 
WRITE (3,700)  EV, TEMP 

700  f  or-mat  ( 1  x 7  ’  ev  suess* '  ,  2 ( el5.  8. 2x )  ,  '  det= '  ,  2  ( el5.  8,  2x )  ) 
f cn*temp 
EVSTORE=EV 

900  formatdx,'  call  to  srid  complete  '7  1fg12.5) 

901  formatdx,'  call  to  bsplcf  complete  ’,lpel2.5) 

902  format(lx7'  call  to  orset  complete  *,lpel2.5) 

903  formatdx,  '  call  to  deset  complete  ',lpel2.5> 

904  format(lx7’  call  to  matrix  complete  ',lpel2.5) 

905  formatdx,'  call  to  depse  complete  ',lpel2.5) 

906  formatdx,'  call  to  sfeval  complete  ',Jpel2.5) 

910  f or ma t ( 1 x , 13el04 2) 

return 


c 

c 

c 


subroutine  matnrm 
IMPLICIT  REAUK8 (A-H, 0-Z) 


PARAMETER (NP  = 
l  n  e  9  p = 

3  nsiF  = 

)  n  d  i  p  = 

5  ntip= 


n  3  d  p  =  1 0 1 , 

nm2p=np-2>  nml F=np-1,  nplp=np+l,  rip2p=np+2» 
neissp=ne««p!|tne^p»  nep=4#ne9F+l  > 

maKPasne*jp*n9dp»  nsap=16#np2p»  rialp=7#ne<;«5,»p> 
ndm2p=ne9P  *np2p  ,  ndm3p=neeis<!tp!|(np2p  , 
n dm23p=3#  nes p  ,  ridm33F=3*neasqp ) 


common/ baksub/aa (ndm2p»  ndm2p) 

c ommon/e3nvc s/e9vs  (np2p>  ne«tp)  ,  usi  (ndm2p)  ,  wa  (ndm2p) 
c ommo n / b d y c d s / b c  1  ( 3 ,  ri e 3 p )  ,bcr(3.neip) 
common/matrix/  a (hp2p»  nesp, ne^F) , 

1  b (np2p, ne^Pi nesp) , 

2  c  (hp2p »  nc-9P,  ne^p) 

c ommon/tnsf rrn/  ar (np2p» ne^Pi neap) ,  ai (hp2p* ne^F,  ne^p)', 

1  br  (hp2p»  neipi  neap)  ,  bi  (np2p.  nesp.  ne«ip)  » 

2  cr(riF-2p>nesp,nesF))  c  i  <np2p,  netFi  ne«ap) » 

3  ars  ( np2p ,  netr  >  netr)  ,  a  i  s  (np2e » rtesp  ,neip)  , 

4  br  s  <np2p,  ne«tP»  ne*tp)  ,  bis  (tip2f  ,  nesp»  neip)  , 

6  crs  ( n  p  2  p  >  neir.  netp)  >cis  (np2p.  ne<ip»  neip) 

common/spvals/phi (maxp) , sevals (nsdFt np2p) 
common/sr ids/nsd»  xlixr»phY9rd(np)»plt9rd(nsdp) 
common/intser /nea, ne^ss,  ndm2» ndm3,  ndm23,  ndm33tmax, itrmax 


ai (np2p» nesp» nesp>% 
bi  (np2p»  ne«tp.  rie«tp)  » 
c  i  ( np2p  >  ne^p ,  nenp )  7 


COMPLEX  ev,evold»e9vs)U'ii)Ciiixe,aaial7phi 
COMPLEX  a, b, c , det, detnrm 
COMPLEX  be  1  *  be r* z norm 


do  10  islindm2 
do  10  j  =  1 1 ndm2 
aa(J7 i)  =  aa < J ,  i  ) #anormi 
10  continue 
c 

return 

c 

entry  setnrm 
anorm  =  ABS(aa(575)) 

PRINT  9000 
WRITE (3, 9000) 

9000  FORMAT (IX, 'ENTER  NORMALIZATION  FACTOR  -  SETNRM') 
READ (5,*)  FACNORM 
WRITE(6,*>  FACNORM 
WRITE (3,*)  FACNORM 
ANORM= ANORM* FACNORM 
anormi  =  1. /anorm 
return 
end 

subroutine  matrix (ev) 

C  IMPLICIT  REAL*8(A-H,0-Z> 


PARAMETER (NP  ■  21, 
J  ne<ip*  1  , 


n  s  i  f  —  1 , 

ridip=  1< 
rit-  i  f=  4, 
n  3  d  f  =  1 0 1  . 

nm2p=np-2,  nml  p=np-l  ,  hf  1  f=hf+1  ,  np2r=n f+2, 
neasaF “negpfcnegp ,  neF=4lnesp+l  > 
nriaxF=Tiegp*nsdP ,  nsap=16#np2p ,  ns  1  p=7»ne«isiF , 
ndm2p=negp*nF  2p ,  ndm3p“neasgp*np2p > 
ndm23p=3)t'nesp  ,  n  dm33p =3  #  n  e  g  s  g  p  ) 


DIMENSION  BTH(NP2P> , PRESS (NP2P) , RHO (NP2P) , 

*  XNU (NP2P) ,VA(NF'2P) ,CS(NP2P) ,  AAAA(NF'2P> ,CCCC(NP2P) , 

*  ALPHA ( NP2P )  ,  00  <  NP2P ) , ALPHP ( NP2P ) , 

*  ALPHR (  NP2P ) , ALPH I ( NP2P ) , ALPHRS ( NP2P ) , ALPH I S ( NP2P ) , 

*  ALPHRP ( NP2P ) , ALPH I P ( NP2P ) 

DIMENSION  FFF (NP2P) ,FFR(NP2P) ,FFI ( NP2P ) ,FFRS(NP2P> , FFIS(NP2P) 

*  FFRP ( NP2P ) , FF I P ( NP2P ) »  FFP (NP2P) 

DIMENSION  GGG(NP2P) , GGSPL <NP2P) ,GGP(NP2P> 
c ommon/bak  sub/aa  < ndm2p *  ndm2p ) 

common/esnvc  s/esvs  (  np2p ,  neip )  1  us  i  ( ndm2p )  ,  wa  <  ndni2p) 
c omiTion/bd vc ds/bc  1  <3?  nesp)  ,  bcr  <3.  n e 3 f ) 


c  ommon/matr ix/ 


c  ommo  n / 1  ns  f r  m/ 


a  (np2p,  nesFine-sp)  , 
b  ( np2p  ,  negp .  near  )  , 
c  <  hp2p , neap 1 neap ) 
ar (np2p. neip. negp) , 
br  (np2p,  negp,  ne«»p)  1 
cr  (np2p»  neap,  nesp)  , 


ai (np2p , negp , nesp ) , 
bi (np2p»  negp, neap) , 
c i <np2p* negp, neap) , 


3  ars (np2p» near, negp) , ais (np2p, negp, negp) , 

)  brs ( hp2p » neap, neip) »  bis (np2p, neip, negp) * 

j  crs (np2p»  neiPi neip) ,  c  is (np2pi neap, negp) 

common/spvals/phi  (tnaxp)  ,  spvals  (nsdf,  np2p> 
c ommon/ grids/ ns d,  xl ,  xr,  physr  d  (np)  1  plt9rd  (nsdp) 
c  ommon/  i  ntser/nM  ,  nesss ,  ndm2,  ndm3,  ndm23,  ndm33,  max  7  i tr  max 

COMPLEX  ev7evoldiesvsiUsi7cmxe7aa7al7Phi 
COMPLEX  a»  b, c  7  det,  detnr m, deval 
COMPLEX  be  1 » be r 7 2 norm 

COMPLEX  AAA, CCC,AAAA,CCCC,FFF, ALPHA, OQ,FFP,ALPHP 


c ommon /nc 1 /nm2, nml , n , n p 1 ,  n p 2 
common/r cl/r (np) , rn (np) 
c  ommon/sac 1 /sa ( nsap ) 

c  ommon/ ph  iOl  /  f-01  (np )  /  phi02/p02  (np )  /  ph  i 03/p03  ( np ) 
common /Phi  1 1 /p 1 1 (np) /phi  12/p  12 (np) /phi 13/p 13 (np ) 
common/p hi21  /f-21  (np)  /phi22/p22 (np)  /  phi23/p23 (np ) 
common/phi il/pil(np)/phii2/pi2(np) 
c  ommon/ phi i3/pi3 (np) /phi i4/pi4 (np) 
c  ommon/ph i4/e ( np2p ) /ph i5/f ( np  2p ) 
common/bndv 1 s/bcOp, bcOp 1 , be  If  1 , be  Ip , be  1 , bcO 
c  ommon/erre 1/err ( np2p ) 
common/serrc  I/amxer  1 ,  amxer-2,  ermax 
common/duml/dl (nr2p) /dum2/d2(np2p) 
c  ommon/c  oef 1 /deva 1 , neva 1 


c 

c#### 

c#### 

c#### 

c#### 


commons  for  spline  integrals 

ssi (nF2, i ) , dsi (nF2, 7, J ) , tsi (np2, 49, k ) 
i “number  of  single  integrals 
j “number  of  double  integrals 
k=number  of  triple  integrals 


( 


c  ommon/ssi  c  1  /  ss  i  (  tip2p  ,  ns  i  f  ) 
c ommon/ dsicl/dsi  (nF-2p,7.ndip) 
c  ommon/ 1  s  i  c  1  /  is  i  <  np2F , 49, nt  3  f ) 

c ommo ri /  nssip/nsii  n s J  (nsip)  ,  nsv  ( 2  >  n  s  i  f-  ) 
common/ndsiF/ndi  ,  ndJ  ( n d i f  )  ,  ndv  (4.  ndip) 
c  ommon  /  n  t  s  i  F-  /  n  t  i  ,  n  t  J  ( n  t  i  f  >  ,  n  t  v  <  6 ,  n  t  i  f-  ) 

equivalence  (sa,sa3) 
dimension  sa3 <4, 4, np2p) 

losical  lintl 

dimension  ibc 1 (nesp) , i bcr (nesp) 

COMPLEX  evsues,ca 
data  1 in  1 1 / . true./ 
data  xlixr/O.il./ 

DATA  GAMMA/ 1. 66667/, ITRMAX/1/ 
data  ibcl, ibcr/nesp*l,nesptl/ 

DATA  AWALL/ 1 . / , XMO/ 1 . 256637E-6/ 

COMPLEX  q,qF,arsq,ssqrt 

the  choices  of  boundary  conditions  are 
=  0  derivative  =  0 
i  b  c 1 , i b  c  r  =  1  function  =  0 
=  2  w.  k.  b. 


if  (lintl)  so  to  1 

•  C*****CALCULATION  OF  COEFFICIENTS 

PHYGRD ( 1 )  =  1 . E-5 

VA ( 1 ) =VA ( 2 ) *FHYGRD ( 1 ) /PHYGRD ( 2 ) 

DO  210  1=1, N 
RR=PHYGRD ( I ) 

CS2=CS(I)**2 

•  VA2=VA(I)**2 
AAA=VA2* (XM/RR) **2-EV 
CCC=CS2*AAA-EV*VA2 
AAAA ( I ) =AAA 

CCCC ( I ) =CCC 

C*  #  *  *  *L AMBDA1 1 /LAMBDA 1 2  CALCULAT I ON 
C  FFF(I)=-RHO(I)/RR* (CCC+2. *EV*CS2) *AAA 

*  /  (EV*EV+CCC#  (XK*XK+  (XM/RR)  **2>  ) 

C*****  (P*)  ' /R  CALCULATION 

GGG ( I ) =RH0 ( I ) *VA2/RR#*2 
C*****CALCULATION  OF  ALPHA 

ALPHA ( I ) =RH0 ( I ) #RR#AAA*CCC/ (EV*EV 
C  *  +CCC* (XK*XK+ (XM/RR) **2> > 

210  CONTINUE 

DO  215  1=1, N 

FFR ( I ) =REAL ( FFF ( I ) ) 

FF I ( I ) = A I MAG ( FFF ( I ) ) 

ALPHR ( I ) =REAL ( ALPHA ( I ) > 
n  ALPH I ( I ) =A I MAG ( ALPH A ( I ) ) 

215  CONTINUE 
C****t 

CALL  DEPSE (FFR ( 1 ) , FFRS ( 1 ) ) 

CALL  REPSP (FFRS ( 1 ) , FFRP ( 1 ) ) 

CALL  DEPSE (FFI (1) ,FFIS(1) ) 

,,  CALL  REPSP (FF IS ( 1 ), FFIP ( 1 ) ) 

C*t**» 


c 

# 

c 

c 


c 

c 

c 

c 


CALL  DEPSE ( GGG ( 1 ) , GGSPL ( 1 ) ) 


CALL  REPSP (GGSPL <  1  )  , GGP < 1 ) ) 

C*»*** 

CALL  DEPSE ( ALPHR ( 1 )  , ALPHRS ( 1 ) ) 

CALL  REPSP (ALPHRS ( 1 ) , ALPHRP < 1 ) ) 

CALL  DEPSE ( ALPH 1(1), ALPH I S ( 1 ) ) 

CALL  REPSP ( ALPH I S ( 1 ) , ALPH I P ( 1 ) ) 

DO  255  1=1, N 

FFP ( I ) =CMPLX (FFRP  < I ) , FFIP(I) >#DELI 
GGP ( I ) =GGP ( I ) *DELI 

255  ALPHP ( I ) =CMPLX (ALPHRP ( I ) , ALPH IP ( I ) ) #DELI 

CttXt* 

DO  220  1=1, N 
RR-PHYGRD  < I ) 

CS2=CS ( I ) **2 
VA2=VA ( I ) **2 
AAA=AAAA ( I ) 

CCC=CCCC ( I ) 

QQ(I)=RR* (FFF(I) *2. /RR* ( 1 . +EV*CS2/CCC> 

*  +RR*GGP ( I ) +RHO ( I) * (4. *VA2* ( 1 . 

*  +EV1CS2/CCC) /F:R**2-AAA)+FFP(I> ) 

220  CONTINUE 

C***«* 

do  2  l = 1 , n 
x  =  f  h  /  9  r d  (  i  ) 


a  ( l ,  1 ,  1 ) 
b  ( i ,  1 ,  1  > 
c (i, 1,  1) 
continue 
return 


ALPHA ( I ) 
ALPHP ( I ) 
OQ  ( I ) 


entry  setbc ( ev ) 


do  90  1 =1 , nei 


i  b  r  n  c  h  =  i  b  c  1  ( i  ) 


so  to  (10,20,30)  ibrnch 
10  continue 

bcl(l,i)  =  (0.,0.) 

b  c 1 ( 2 , i )  =  (l.,0.) 

®  b  c 1 ( 3 , i )  =  (0.,0.) 

so  to  90 
20  continue 

be  1  ( 1 ,  i  )  =  <1.,0. ) 

be  1 (2,  i  )  =  (0.,0.) 

be  1 (3, i )  =  <0.,0.) 

9  o  t  o  90 
30  continue 

s  =  c(l,i,i)*bcO 

«tp  =  (c  <1,  i,  i)*bcOF+c  (2,  i,  i)*t.cOFl)*deIi 


a  r  9  s  =  9 

if  (REAL  (ar-99)  .  le.O.  .  and .  A I  MAG  (ars9)  .  gt.O.  ) 
<•  959rt  =  S0RT(ar99> 

if (REAL (9) . lt.O. .and. A I MAG ( 9 ) . e9 . 0. )  9E9rt 
be  1 < 1 , i )  =-. 25*9f/9+ (0.,l.)*959rt 
bcl (2, i)  =  DCMPLX (1. ,0. ) 
b c 1 ( 3 , i )  *  (0.,0.) 


ar?R=conJ 9 (ar  99 ) 

=  -  9  s  9  rt 


90  continue 

do  92  i.  =  l  ,  ne9 


i  t'  r  n  c  h  =  i  b  c  r  (  i  )  +  1 

so  to  (12,22,32)  ibrnch 
continue 
b  c  r  ( 1 ,  i  ) 
b  c  r  ( 2 ,  i ) 
b  c  r  <  3 ,  i  ) 
so  to  92 
continue 
b  c  r  ( 1 ,  i  ) 
b  c  r  ( 2 ,  i  ) 
bcr (3, i ) 
so  to  92 
continue 

s  =  c  ( n  p  2 ,  i ,  i  )  *  b  c  1 

se  =  (c(nF2,i,i)*bclF+c(npl,i,i)*bclFl)*deli 
a  r-  3  s  =  9 

i  f  (REAL  (ars9).le.0..and.  A I  MAG  ( a  r  9  s )  .  g  t .  0.  >  ar99=conj9(ar99) 
sssrt  =  SQRT(arss) 

i f  (REAL < s ). 1 1 . 0. . and . AI MAG ( s >. es . 0. )  sssrt  =  -sssrt 
be r  < 1 , i )  =-. 25*sf/s- <0. , 1. ) *sssrt 
b c  r ( 2 , i )  =  DCMPLX (1.  ,0. ) 
b  c  r  ( 3 ,  i  )  =  ( 0 .  ,  0 .  ) 

continue 
return 
continue 
WRITE (6, 1000) 

WRITE (3, 1000) 

read  <5, * )  ev9ues,deval,neval,xk,xm, gamma , ibc 1 , iber, x I , xr , 
i trmax 

wr i te (6, * )  evsues , deva 1 , neva 1 , xW , xm,  gamma  ,ibcl,ibcr,xl,xr, 
i  trmax 

wr i t e (3, * )  evsues , deva 1 , neva 1 , xk , xm, gamma , ibcl, ibcr,xl, xr» 
i trmax 

ev  =  evsues 
REWIND  30 

READ ( 30 )  AW ALL , BTH , RHO , PRESS , VA , CS 
C * * * * * NORMAL I Z AT I ON 
VAA=VA(N) 

RHQN=RHO ( 1 ) 

WRITE (6, 1001) 

WRITE (3, 1001) 

DO  800  1=1, N 

VA ( I ) =VA ( I ) /VAA 

CS ( I ) =C5 ( I ) /VAA 

RHO ( I ) =RHO ( I ) /RHON 

WRITE (6,  1002) I , VA ( I ) , CS ( 1 ) , RHO ( I ) 

WRITE (3, 1002) I , VA ( I ) , CS ( I ) , RHO ( I ) 

800  CONTINUE 
C**##* 

deli  =  1 . / ( xr-x 1 ) 

1 i nt 1  =  . f a  1 se. 


= 

(0. ,0. ) 

= 

( 1 .  ,  0.  ) 

= 

(0. ,0. ) 

II 

(1. ,0. ) 

II 

(0. ,0. ) 

II 

(0. ,0. ) 

12 


92 

1 


1000  format ( lx, ' enter  namelist  data'/ 

*  3X ,  ' EVGUES , DEV AL , NE VAL , XK, XM, GAMMA, IBCL, I  BCR, XL,  XR,  ’ , 

*  ' i trmax ' / 

*  3X , ' EVGUES=( OMEGA* AWALL/ VAA) **2, XK=K*AWALL ' ) 

1001  FORMAT (3X, 'NORMALIZED  EQUILIBRIUM  PARAMETERS'/ 

*  8X,  ' I' , 10X,  'VA(I)  ’ , 10X,  'CS(I)  ' ,9X,  ' RHO ( I )  ' ) 

1002  FORMAT  <3X, I5.3E15.5) 
return 


subroutine  nrmfcn 
IMPLICIT  REAL #8 (A— Hi Q— Z ) 


-ei senf unc t i ons  normalized  such  th3t  the  value  of  the 
central  ei senf unc t i on  is  (l.,0.)  at  x=0. 


PARAMETER (NP  = 
’  n  e  a  p  = 


n  s 1 p  =  li 
ndip=  It 
ntip=  4, 
ri3dp=101  > 

nm2p=np-2->  nmlp=np-l ,  nplp=np+l  ,  np2p=np+2, 
neasap=neap*neap  >  nep=4#neap+l . 
maxp=neap#nsdp, nsap=16*np2p, nalp=7*neasap, 
n d m 2 f- = n e a f- # n p 2 p ,  n d m 3 f  =  n e a s a f •* n p 2 f- , 
ndm23p-3#neap»  ndm33p=3#neasap  > 


c ommon/bak sub/aa (ndm2p» ndm2p) 

common/esnvcs/esvs (np2p,  neap) >usi (ndm2p) , wa (ndm2p> 
comtrion/bdvcds/bcl  (3,  neip)  ,  bcr  <3.  ne«ip) 
common /matrix/  a (np2p, neap, neap) , 

1  . b <np2pi nespt neip) , 

2  c  (np2p,  ne^Pi  neap) 

comman/tnsf rm/  ar  (op2p,  neiF,  neap)  ,  ai  (rip2fi  ne^F-,  neap)  , 

1  br(np2p»nesp»nesp)»  bi  (np-2p»  neap,  neap)  , 

2  cr <np2p» neap, neap) ,  c i (hp2p, neap, neap) , 

3  ars (np2p, neap, neap) , ais (np2p, neap, neap) , 

4  brs  <nr2p, neap, neap) , bis  <np2p*  neap, neap) , 

6  crs (np2p *  neap, neap) , c is (np2p , neap,  neap) 

common/spvals/phi (maxp) , SF-vals  <n9dp,  np2p) 
common/srids/nsd, xl , xr  , phvsr-d (np) ,  pltsr-d (nsdp) 
common/intser/nea, neasa, ndm2, ndm3, ndm23, ndm33, max, itrmax 
c  o  m  m  o  n  /  n  o  r  m  /  z  n  o  r  m 


COMPLEX  ev,evold,esvs,usi,cmxe,aa,al,phi 
COMPLEX  a, b, c , det, detnrm 
COMPLEX  be  1 , bcr , z norm 


c  ommo n / n  c 1 / nm2 , nm 1 , n , n  f 1 , n  p  2 
common/rc 1/r  (np)  ,  rn  (np) 
common/sacl/sa (nsap) 

common/phiOl /pOI (np) /phi 02/p02 (np) /phi03/p03 (np) 
common/phi Il/pll(np)/phil2/pl2(np)/phil3/pl3(np) 
common/Phi21/p21 ( np ) /phi22/p22 (np) /phi23/p23 (np  > 
c  ommo n/ phi i 1/pil (np)/phii2/pi2(np) 
c  ommon/phi i3/pi3(np)/phii4/pi4(np) 
common/phi4/e ( hf2p ) /phi5/f (np2p> 
common/bndvl s/bcOp, be Op  1 , be 1f1 , be  1 p,  be  1 , bcO 
common/errc 1/err (np2p) 
common/ ser-r-c  1  /amxer  1 ,  amxer-2,  er-max 
common/ duml/dl (np2p ) /dum2/d2 (op2p) 


c 

c 

c#### 

c#### 

c#### 

c#### 


c  ommons 


spline  intesrals 


ss  i  ( np2,  i  )  ,  ds  i  (np2,  7,  J  ) ,  ts  i  ( of  2, 49, 1: ) 
i=n umber  of  sinsle  intesrals 


.i -number 
k*=number 


double 

triple 


i  ntesr-a  1  s 
i  ntesr-a  1  s 


c  ommon/ ss  i  c  1  /  ss  i  (  np2p  7  ns  i  p  ) 
c  ommon /dsicl/dsi  <  n f  2f-  ,  7 ,  n d  i  p  ) 
c  ommon/ 1 s i c 1 / ts i ( of  2p . 49. n t i f ) 


c  omiiion/nss i  p /ns i  >  nsJ  ( ns i f ) . ns v (2. nsi p ) 
c ommon/ nds i p/nd i >  ndJ  ( rid  i  f-  )  ,  nd  v  < 4 .  nd  i  p  ) 
c ommon/nts i p/nt i , ntJ (ntip) . ntv  <6. ntip) 

equivalence  (sa.sa3) 
dimension  sa3 <4 . 4 . np2p ) 


do  60  i ter=l . itrmax 

call  leqtlc (aa.  ndm2? ndm2. usi 7 1 . ndm272. us. det. ier) 


do  50  i=l7neq 
do  50  J  =  1 . n p 2 

esvsfj. i)  =  us i ( i+ ( J -1 ) Inei ) 
50  continue 

-  diagnostic  print  - 

do  40  i=l.nes 
do  40  J=l>np2 
WRITE (07 100)  e s  v s  <  J ,  i  ) 

40  continue 
60  continue 


100  f ormat ( 1 x 7 2el2. 5) 

101  format (lx, 13el0. 2) 


-  setup  eigenfunctions  - 

2  continue 

do  32  i  =  1 7  ma x 
phi ( i )  =  ( 0 . 7  0 . ) 

32  continue 

do  22  l!=l7nei 
do  22  i=l.ngd 
index*5  (k-i>  #nsd+i 
do  22  J=l7np2 

phi (index)=phi (index)+esvs(J7k)*spvals(i7J) 
22  continue 

normalize  eigenfunctions 

znorm=phi ( 1 ) 
do  24  i=2. max 

if (ABS ( phi ( i ) ) . st . ABS <z norm) )  znorm-ph i ( i ) 
24  continue 

IF (ABS  <  ZNORM) . LT . 1 . E-20)  ZN0RM=(1. ,0. > 
znorm=cmplx (abs (znorm) 7  0.  ) 


return 


c 


subroutine-  of- set 
IMPLICIT  REAL*8(A-H,0-Z) 


c 

c 


PARAMETER (NP  =21, 


o 

n  e  q  p  = 

1 , 

3 

ns  i  f-= 

1, 

4 

n  d  i  F-  = 

1 , 

5 

nt  i  p= 

4, 

6 

n  s  d  F-  = 

101 , 

-7 

/ 

nm2p= 

np-2 ,  nmlp=np-l ,  riF- 1  p=mp+1  ,  riP2p=np+2, 

7 

neqsqp=neqp$neqp , nep=4*neqp+l , 

3 

maxp= 

neqp#nsdp,  nsap=16#np2p,  na  1  p=7#neqsqp 

4 

ndm2p 

=neqp#np2p, ndm3p=neqsqp*np2p , 

5 

ndm23p=3#neqp,  ndm33p=3Htneqsqp ) 

c 

c 

c  o mmo  n  /  n c  1  /  nm2 ,  nm  1 ,  n ,  n p  1 ,  n  f-2 
common/rc  1/r  (riF-)  ,  rn  (np) 
c  ommon/sac 1 /sa  ( nsap ) 

c  ommon/phiOl/pOl  <  hp )  /phi02/F-02  (rip)  /ph i03/p03  (np ) 
common/phi  1 1  /  f-  1 1  (rip )  /  f  hi  12/p  12  (of  )  /  phi  13/p  13 (of  ) 
comrr.on/phi21/F-21  (np)  /phi22/p22(np)  /phi23/p23 (np) 
common/phi  i  1/pi  1  (np)  /f  hi  i2/pi2  (rip) 
common/phi i3/pi3 (np ) /phi i4/pi4 (np ) 
common/phi4/e(riF-2p>  /phi5/f  (np2p) 
c  ommon/bndvl s/bcOp,  be  Op  1 ,  be  1  f-  1 ,  be  Ip  ,  be  1 ,  bcO 
common/er-rcl/err  (np2p) 
common/serrc 1 /amxerl ,  amxer2»  er-max 
common/duml /dl  (nF-2p)  /d um2/d2(np2e) 
c 

c  commons  for  spline  intesrals 

c 

c####  ssi  ( np2,  i  )  ,  dsi  ( riF-2, 7,  J  )  ,  ts  i  (np2,  49,  k  ) 
c####  i=number  of  sinsle  intesrals 

c####  J  “number  of  double  intesrals. 

c####  k=number-  of  triple  intesrals 

c 

c  ommo  n  /  s  s  i  c  1  /  s  s  i  ( n p2f-  ,  n  s  i  f-  ) 
c ommon/ds  i  c  1  /d.s  i  ( np2p ,  7,  nd ip ) 
common/tsic 1/tsi (np2p , 49, nt i p> 


common /ns si  p/ns i , nsJ (ns  ip) , nsv  <2, ns  ip) 
common/ndsip/ndi , ndJ (ndip) , ndv (4, ndip) 
c  ommon/nt s i p/nt i , nt J ( nt i p ) »ntv(6,ntip) 


c 

v 


-  t  /. 


c* 

\%\ 

V- 

a  - 

i- 


equivalence  (sa,sa3) 
dimension  sa3 (4, 4, ne2p) 


data  nsJ/nsip#0/ 
data  ndJ /ndip*0/ 
data  ntJ/ritiP*0/ 
data  nsv/0,0/ 
data  ndv/ 1 , 0, 1 , 0/ 
data  ntv/O, 0, 1 , 0, 1 , 0, 
2  1,0, 1,0, 0,0, 

3  0,0, 1,0, 0,0, 

4  0,0, 0,0, 0,0/ 
return 

end 

subroutirie  Fltfen 


IMPLICIT  REAL*a<A-H,0-Z) 


PARAMETER <NP  =21, 


neif'=  1  , 
n  5  i  f  =  1 , 

n  d  i  f  =  1 , 

ntip=  4, 
n  9  d  p  =  1 0 1 , 

nm2p=np-2,  nmlp=np-l ,  nplp=np+l ,  riP2p=np+2, 
neps9p=ne9p#nepp ,  nep=4#nepr+l , 
maxp=ne9F  # ns d p , n sa r= 16# np2p ,  na 1  p=7# n e  p s q f  , 
n  dm2p=n  e  9 f-  # n f-2f-  ,  n  dm3 p =n  e p  s  q  p # n  f-2f  , 
ndm23p=3#neqp,  ndm33p=3#nepsqp  ) 


c  ommon/bak sub/aa  ( ndm2p ,  ndm2p ) 

c ommon/esrivc s/b3 vs  (riF 2f  ,  riB9F  >  ,  us i  < ndm2p )  ,  wa  (ndm2p) 
common/bdYcds/bc  1  (3,  nepp)  ,  bcr  <3,  nepp) 
common/matrix/  a ( np2p , nepp , nepp) , 

1  b (np2p, nepp, nepp) , 

2  c (np2p» nepp, nepp) 

common/tnsf rm/  ar-  (np2p,  nepp,  nepp)  ,  ai  ( hp2p ,  nepp,  nepp)  , 

1  br (np2p, nepp, nesp) »  bi (np2e» nepp? nepp) , 

2  cr  (hp2p» nepp, nepp) ,  ci (dp2p, nepp, nepp) , 

3  ars (np2p, nepp, nepp) , ais (np2p, nepp, nepp) , 

4  brs (np2p, nepp, nepp) , bis  <np2p>  nepp, nepp) , 

b  crs (np2p, nepp, nepp) , c is (np2p,  nePF  ,  nepp) 

common/sevals/phi  (maxp)  ,  spvals  (nsdp,  np2p) 
common/srids/nsd, xl , xr, phvsrd (np) , pltsrd (nsdp) 
common/intser/nep,  nepsp,  ndm2,  ndm3,  ndm23,  ndm33,max,  itrmax 
common/ norm/s norm 

COMPLEX  ev,evold,esvs,usi,cmxe,aa,al,F-hi 
COMPLEX  a , b, c , det , detnrm 
COMPLEX  be  1 , be r , znorm 


dimension  f  (nsdp)  ,  F-ltr  (nsdp)  ,  F-lti  (nsdp) 

common/nc  1  /nm2,  nml ,  n»  nF-1 ,  hp2 

COMPLEX  temp ( np2p , nepp, nepp) 

REWIND  15 
do  2  i=l,nsd 

2  Pi tsrd ( i ) = (xr-xl ) #p1 tsrd ( i ) +xl 
do  100  m= 1 , 4 
DO  3  1=1, NOD 
PHI (I)  =  (0.  ,0.  ) 
so  t o ( 1 0 , 20 , 30 , 40 ) m 

10  continue 

do  11  i=l,nep 
do  11  J=l,nep 
do  11  k=l,np2 
t  emp ( k , J ,  i )  =  a ( k , J , i ) 

11  continue 
so  to  50 

20  continue 

do  21  i=l,nep 
do  21  J=l,nep 
do  21  k=l,np2 
temp (k, j , i)  =  b  (k , J , i) 

21  continue 


so  to  50 
30  continue 

do  31  i=l,nes 


do  31  J=l,nes 
do  31  k  =  1 ,  n f- 2 
t  erne  ( k  ,  j  ,  i  )  =  c  ( k  »  j  »  i  ) 

31  continue 
so  to  50 

40  continue 

do  41  J=l,nes 
do  41  k=l,np2 

temp  ( k  ,  1 ,  J  )  =  esvs  (k  » J  )  /znornt 

41  continue 
50  continue 

do  60  i=l,nes 

do  61  J=l,nes 
do  62  k  =  1 ,  n  3  d 
do  62  1=1, ne2 

phi(k)  =  phi (k)+temp (1, i» J ) *spvals(k» 1) 

62  continue 

do  64  mm=l,nsd 

pltr (mm)  =  REAL ( Ph i (mm) ) 

plti(mm)  =  AIMAG (phi (mm) ) 

64  continue 

call  extr-ma  (Pltr,  pi  ti ,  n9d»  vmin,  vmax) 
xstp  =  .  1  *  (xr-x 1 ) 
vstp  =  . 2* ( vmax-vmin) 

WR I TE ( 90 )  NGD , PLTGRD , PLTR , PLT I , YM I N , YMAX 
call  title( ' prosram  esv 1 et$ ' , -100, ' doma i n ’ > 6, 

2  ' amp  1 i tude ' , 9, 9. , 6. ) 

call  9raf(pltsrd(l),xstp,pltsrd(nsd),Ymin,Yr.tP,Ymax) 
call  curve( p Its rd,  pltr  ,  ri3d,0) 
call  curve(Pltsrd,Plti,ri3d,0) 
call  frame 

do  68  nn=i,n3d 
68  Phi (nn)  =  (0.  ,  0. ) 

61  continue 

if(m.es.4)  9oto  101 
60  continue 
100  continue 


continue 

return 

end 

subroutine  externa  ( 1 1  *  t2->  n>xl.x2) 

IMPLICIT  REAL#8(A-H,0-Z) 

dimension  tl(i),t2(l) 

xl=tl ( 1 ) 

x2=t 1  (1) 

do  1  i  =  l ,  n 

if (xl. st. tl  (i  ) )  xl=t  1 ( i ) 

if (xl. st.t2(i) )  x 1 =t2 ( i ) 

if <x2. It. tl  (i) )  x2=t 1 ( i ) 

if (x2. It. t2(i) )  x2=t2 ( i ) 

continue 

return 

end 

prosr  am  s Ysode 
IMPLICIT  REAL*8(A-H,0-Z) 


PARAMETER <NP  =21, 
neap=  1, 
nsip=  1, 
ndip=  1, 
ntiF=  4, 
nsdp=101, 

nm2p=np-2,  nmlp=np-l ,  nr lF=np+l ,  nr2r=nr+2, 
neasap=neaF#neap,  nep=4*neap+l , 
maxF=neap*nsdP,  nsap=16*np2r »  na lp=7#neasaF, 
ndm2p=neap*np2p,  ndm3p=neasap*np2p > 
i  ndm23p=3#ne*»p,  ndm33p=3#neasap) 

common  / 1  osi  c  1  / 1  intal ,  1  dtnrm,  lrerdr 
1 03 i ca  1  1 inta 1 , 1 dtnrm, lrerdr 

common /baksub/aa <ndm2r, ndm2p) 

common/ e9nvc  s/esvs (nr2r , nesr ) , usi (ndm2r  > , wa ( ndm2p) 
common/bdYcds/bc 1 (3, neap) , bcr (3, nesr) 
common/matrix/  a ( nr2r , neap,  neip) , 

b  <np2p»  near, nesr) , 

!  c <np2p> neap, neap) 

c ommon/tnsf rm/  ar (np2r, neap, neap) ,  ai (nr2p*  near, neap) , 
br  <Yip2p ,  neap,  near)  ,  bi  (nr2p,  neap,  neap)  , 

!  cr (np2p, near, near) ,  c i (np2p, near, near) , 

!  ars (np2r, neap, neap) , a  is (nr2p, neap, neap) , 

brs (np2p»  neap, neap) , bis (np2r, near, near) , 
i  cr s  <np2p, neap, neap ), c is (nr2p, neap , neap) 

common/spvals/phi (maxr) , srvals (ns dp, np2p) 
common/srids/nsd, xl,xr,phYsrd(np),pltsrd(nsdp) 
common/intser/nea, neasa, ndm2, ndm3, ndm23, ndm33,max, itrrnax 

COMPLEX  ev,evold,esvs,usi,cmxe,aa,al,Phi, ELF 
COMPLEX  a, b, c , det, detnrm, deval 
COMPLEX  bcl,  bcr,. z norm 


common/ nc 1 /nm2, nml , n, nr  1 , nr  2 
c ommon/rc 1 /r ( np) , rn  (nr) 
c  ommon/sac 1 /sa (nsap ) 

common/phiOl/pOl (np) /phi02/p02(np) /phi03/r03(nr) 
c ommon/Phi  1 1  /f-  1 1  (np )  /ph i  12/p  12  ( np )  /  phi  13/p  13  (np ) 
c  ommo n/ p h i 2 1 / f  2 1 ( n r ) / f  b i 22/ r 22  <  n r ) / p h i 23/ f  23 ( n p ) 


u  a  u  u  u 


c#### 

c#### 

c#### 

c#### 


c  ornmo n / f  h 1 1 1 / p 1 1 <  n  f ) / p  h 1 1 2/ p 1 2  ( n p  > 
c  ommon/ph i  i3/p  i3  < riF  )  /f  h  i  i4/p  i4  ( np ) 
c  ommon/ph i 4/ e ( np2p  > /ph i5/ f ( hp2p ) 
c  ommon / bnd  v 1 s / b  cOp  ?  be  Op  1 ?  b  c 1 p 1 ?  be  1 f  ?  be  1 ?  bcO 
c ornmon/er  r  c  1  /err  ( np2p ) 
c  ornmon/ser  r  c  1  /aitixer  1 ?  amxer  2 ?  er-rna x 
c  ornmo n /d  urn  1  /d  1  (of  2f-  )  /dum2/d2  (np2p) 
common/c oef 1 /deva 1 ?  neva 1 

commons  for  spline  integrals 

ssi  (of  2?  i  )  >dsi  (  np2?  7.  J  )  ?  t  s  i  (np2?  49?  k  ) 
i=numt<er  of  sin9le  integrals 
J=number  of  double  integrals 
k=number-  of  triple  integrals 

common/ssicl/ssi (np2p?  nsip) 
common/dsic 1/dsi (np2p,  7.  ndie) 
common/ tsic 1/tsi (hf2p ? 49?  ntip ) 

common/nssip/nsi ?  nsj (nsip) ? nsv  <2?  nsip) 
common/nds i p/ndi ?  ndJ (ndip) ?ndv(4?ndip) 
common/ ntsip/nti » ntj (ntip) ?  ntv  <6, ntip) 
COMMON/ESTORE/EVSTORE 

equivalence  (sa>sa3) 
dimension  sa3 <4, 4. np2e) 


data  mri2)  nml ,  n>  riPl ?  np2/nm2p»  nmlp?  hp»  ppIp,  np2e/ 
data  nsi?ndi?nti/nsie?ndip?ntip/ 

data  nes,  nad?  neqss?  ne<max,  nsa/neiF  .  nsdp»  neqsqp?  nep?maxp?  nsap/ 
data  nal  7  ndm27  ndm37  ndm237  ndm33/nalP7  ndm2p?  ndm3p,  ndrn23p.  ndm33p/ 

data  Iintal7ldtnrm7lrerdr/.true.7.false.7.false./ 


COMPLEX  fstdet, fcn?EVSTORE 
external  fen 

call  1 ink ( ' uni t3= (output 7  create?  text ) ?  prints?  uni t59=ttv// ’ ) 

CALL  TEKALL <4014, 120,0? 0?0) 

CALL  BGNPL(O) 

print  100 
WRITE (3, 100) 

f ormat ( 1 x? ' f irst  call  to  fen  initiated') 
f stdet=fcn (ev) 

EV=EVSTORE 
print  101 
WRITE  <3? 101) 

f ormat ( lx? ' first  call  to  fen  complete') 

if (neva 1 . ne. 0)  soto  200 
print  102 
WRITE (3, 102) 

f ormat < 1 x ? ' ca 1 1  to  cevalf  initiated') 

ELF=  cevalf (ev ? fen) 
print  103 
WRITE (3? 103) 

f ormat ( 1 x ? ' ca 1 1  to  cevalf  complete’) 


call  nrmfcn 


call  Fltfcri 

9oto  300 
200  continue 

call  endrl  (0) 
call  donee  1 
rewind  90 
write(90>  neval 
do  10  nn=l i neva 1 
det=fcn (ev) 
write (90)  ev,det 
CALL  NRMFCN 
call  plt.fcn 
ev=ev+deva 1 
continue 
endfile  90 
300  continue 
stop  999 
end 

subroutine  sprnt  (ar in, irouti imax) 
IMPLICIT  REAL*8(A-H,0-Z) 


input 

arin:  beginning  location  of  floating  point  numbers  to  process 

irouts  beginning  location  to  store  compact  bed  representation 
of  floating  point  numbers 

imax:  total  no.  of  successive  floating  point  numbers  to  process 

purpose 

to  convert  n=imax  successive  floating  point  numbers  into  a  compact 
bed  representation  in  the  e-format  and  store  them  in  n  successive 
teiriForarv  locations  beginning  at  irout  which  are  later  to  be 
output  in  r-6  format. 

for  example  the  floating  point  number  -1.7465e+25  is  printed 
as  -175+g  where  the  decimal  point  is  assumed  between  the  minus 
sisn  and  the  di9it  1.  the  2  digit  exponent  is  converted  to  a  single 
character  bv  a  table  lookup  on  (0, 1  * . . . » 9. a* b, . . . 2 ) ® <0, 1 , 2, . . . 35) 

if  an  i.  or  r.  or  #.  is  returned  by  subroutine  etYPe.  it  is  printed 

if  the  exponent  is  greater  than  +35  or  less  than  -35  a  *  is 
printed 


REAL* 16  aout 


dimension  arin( 

2 )  ,  i  r  o  u  t  ( 2 )  , 

iex (73) 

data  iex/4H  -* 

,4H  - 

* 

X 

<4“ 

N 

y,4H  -x 

,4H  - 

i 

X 

3 

V 

,4H 

-u  ,  4H 

-t , 4H 

—  5 

, 4H 

-r  i  4H 

-s  ,  4H 

-r,4H 

i 

o 

I 

-n > 4H 

1 

3 

w 

X 

-1 

, 4H 

-k  >  4H 

-J  ,4H 

-i,4H 

-h 

,4H 

— g  i  4H 

-f  ,4H 

-e,4H 

-d » 4H 

-c,4H 

X 

<4- 

jUj 

\ 

-a 

,4H 

-9*  4H 

-8, 

4H 

-7, 4H  - 

6,  4H 

—5,  4H 

-4,4H  - 

3,  4H 

4- 

» 4H 

-1,4H 

+0.4H 

+ 1  *  4H 

+2 

,4H 

+3*  4H 

+4, 4H 

+5,  4H 

+6,4H 

X 

* 

ir 

+8>  4H 

+9 

nnnnnnnnnnnnn 


.  4H 

+a ,  4H 

+  b .  4H 

+c,4H 

+  d ,  4H 

+e 

,  4H 

+  f  ,4H 

+  S.4H 

+  h ,  4H 

+  i,4H 

+J  7  4H 

+  k  ,  4H 

+  1 

,  4H 

+m » 4H 

+  n .  4H 

+  o .  4H 

+p 

.  4H 

+  s.4H 

+  r,4H 

+  s .  4H 

X 

r» 

y 

+ 

+u,  4H 

+  v,4H 

+w 

,  4H 

+  x,  4H 

+  v .  4H 

+z ,  4H 

+  #/ 

do  100 

i  =  1 . 

imax 

convert  arin(i)  to  e9.2  format 
call  zcetoa (aoutiO.arin(i) i9i2) 
find  exponent 

call  zciatob(aout>&i8iind’2) 
check  if  index  within  bounds. 

ind-maxO  (minO  ( i  n d  ■>  35 )  ,  -37) 
pick  up  exponent. 

irout(i)  =  iex(ind+38) 

pick  up  si9n  bit  and  first  significant  disit. 

call  zmovechr-  (  irout  (  i  )  ,  0.  aout .  0. 2) 
pick  up  second  and  third  significant  digits, 
call  z  mo  v  ec  h  r ( i r  o  u  t ( i )  .2. aout. 3. 2) 

Cl 00  continue 
return 
end 

subroutine  prnt  (arraYi idim.Jdim. imin. imax. Jmin. Jmax) 

C  IMPLICIT  REAL*8(A-H,a-Z) 

PARAMETER  (icr  =  19) 
c .. programmer s  J  breazeal 
c  . .  date:  1/10/75 

c..  things  to  consider  in  using  this  subroutine 

•  c  set  value  of  icr 

c  set  value  of  nout 

c  add  declarative  for  1cm  if  array  is  in  1cm 
c  add  declarative  value  i d im» j dim* imin» imax > Jmi n» Jmax 
c  this  routine  calls  srrnt  (J  breazeal)  and  orderlib  routines 

•  c..this  routine  is  fashioned  after  the  'nrl'  prnt  routine  (see  d  ander-so 

c..this  routine  will  print,  the  2-d  array  specified.  a  partial 

c.. printout  of  the  array  is  obtained  by  setting  imin. imax.  and 

c .. Jmi n» Jmax.  the  maximum  points  per  line  is  19  in  the  i  direction. 

c.. there  is  no  limit  in  the  J  direction.  if  more  than  19  points 

c..span  the  printout  in  the  i  direction,  the  printout  occurs  in 

c.. partial  blocks,  where  i  =  imin  to  imin+19,for  Jmin  to  Jmax* 

c..  and  next:  i=imin+19  to  imin+39,for  Jmin  to  Jmax!  etc. 

c . . usuase: 

c  call  prnt  (array ( 1 . 1 )» idim. J dim. imin. imax. Jmin. Jmax) 


c.. array:  array  to  be  printed 
c..idim:  dimension  in  i 
c..Jdim:  dimension  in  J 

/r  c ..  imin.  imax:  low  and  high  range  of  index  i 
c .. Jmin. Jmax:  low  and  hieh  range  of  index  J 
c..icr!  max  no.  of  points  to  print  per  line 


PARAMETER  (nout  =  3) 

dimension  array  ( idim. J dim) » sparay ( icr ) 


n  n 


do  25  it'll:  =  imi n , imax , i c r 

i2  =  i2  +  i  c  r  +  imin  -  1 
if  ( i2. st. imsx)  i2  =  imax 
in  urn  =  i2  -  iblk  +  1 
write  (noutiSl)  <(i,i=iblk,i2)) 

do  20  J  =  J  m  i  n ,  j ma x 


call  sprnt  (arrav  (iblhi  J  )  i  SF-aravi  inurri) 
write  (nout,50)  (J  ,  (spar-av  (  i  )  ,  i  =  l ,  in  urn)  ) 
20  continue 

25  continue 


50  format  <3h  J  = , i  2 ,  19a6) 

51  format  ('  i  ='  , 19i6) 
return 

end 

subroutine  legtlc  (a»  n»  ia,  b».m,  ib,  iJob,  wa>  det,  ier) 
IMPLICIT  REAL#8 ( A-H ,  0— Z ) 


c 

c-lestlc - 

c 

c  function 

c 

c 

c  usage 

c  PARAMETERS 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

cm 

c 

c 

c 

c 


s/d 


a 


n 

ia 


b 


m 

ib 

i  J  ob 


wa 

det 
i  er 


library  1 - 

-  matrix  decomposition,  linear  equation 

solution  -  space  economiser  solution  - 
COMPLEX  matrices 

-  call  leatlc  (a,  ri,  ia,  b»  m,  ib ,  i  J  ob,  wa ,  det ,  i  er) 

-  input  COMPLEX  matrix  of  dimension  n  bv  n 

containing  the  COMPLEX  coefficients  of  the 
equation  ax  =  b. 

on  output,  a  contains  the  1-u  decomposition  of 
a  rowwise  permutation  of  the  input  matrix  a. 

-  order  of  matrix  a.  (input) 

-  row  dimension  of  a  as  specified  in  the  calling 

program,  (input) 

-  input  COMPLEX  matrix  of  dimension  n  by  m 

containing  the  m  COMPLEX  valued  right  hand 
sides  of  the  equation  ax  =  b. 
on  output,  the  solution  matrix  x  replaces  b. 
if  i  J  ob=l ,  b  is  not  used. 

-  number  of  right  hand  sides  (columns  in  b). 

(input) 

-  row  dimension  of  b  as  specified  in  the 

calling  program,  (input) 

-  input  option  PARAMETER.  iJob=i  implies  when 

i=0,  factor  the  matrix  and  solve  the 
eaua t i on  ax=b. 
i = 1 ,  factor  the  matrix  a. 
i=2,  solve  the  equation  ax=b.  this 

option  implies  that  legtlc  has  already 
been  called  using  iJob=0  or  1  so  that 
the  matrix  has  already  been  factored.  in 
this  case  output  matrix  a  must  have  been 
saved  for  reuse  in  the  call  to  lestlc. 

-  work  area  of  length  n  containing  the  Pivot 

indices. 

-  determinant  of  matrix  a. 

-  error  PARAMETER, 
terminal  error-  *  128+n 

n=l  indicates  that  matrix  a  is 

algorithmically  singular.  (see  the 


*.  v 


F-recisi  on 

reqd.  ims  1 '  r  out  i r.es 
language 


chapter  ].  prelude) 
s i n3l e/doub 1 e 
uertst 
fortran 


I  c 

cm 

cm 

cm 

c 


latest  revision  -  January  8,1973 

addition  of  determinant  calculation 


mav  17,1979 


subroutine  leqtlc  (a , n , ia , b, m» i b, i J ob, wa , det , i er ) 
IMPLICIT  REAL #8 ( A-H, Q-Z ) 


dimension 
double  precision 


a  (ia,  1) ,  b  (it.,  1) ,  wa  (1)  ,  t  (2) 
p , q  >  z  er  o , one, wa , t , r  n » b i 9 


COMPLEX 

COMPLEX 

equivalence 

data 

data 


a, b, sum, temp 

a , b , sum, temp , det 
( sum, t ( 1 ) ) 

zero/O. dO/,on e/1. dO/ 
zero/0.0/, on e/1.0/ 


initialization 


ier  =  0 
if  (iJob  .  eq. 
r  n  =  n 


!)  9  0  to  75 


find  equilibration  factor 


do  10  i = 1 , n 
bis  =  zero 
do  5  J  =  1 , n 

temp  =  a(i»J) 
p  =  cdabs(ternp) 
p  =  ABS(temF) 
if  (p  . 9 1 •  bi9)  bi9  =  p 


det  =  (1. ,0. ) 
ssn  =  1. 

5  continue 

if  (bis  .eq.  zero)  90  to  105 
wa  ( i )  =  one/bis  . 

10  continue 

1-u  decomposition 

do  70  J  =  l,n  . 
j  m  1  =  j  —  1 

if  (Jml  .It.  1)  90  to  25 

compute  u(i,J>,  i=l,...,J 

do  20  i=l,Jml 
sum  =  a  ( i , J ) 
iml  =  i - 1 

if  (iml  .It.  1)  so  to  20 
do  15  k=l , iml 

sum  =  sum-a ( i , k > #a ( k , J ) 

15  continue 

a(i,J)  =  sum 
20  continue 

25  p  =  zero 

compute  u ( J , J )  and  l(i,J) 

do  45  i=J , n 

sum  =  a  (  i , J ) 

if  (Jml  .It.  1)  90  to  40 
'  do  35  k = 1 , J m 1 

sum  =  sum-a ( i , k ) *a ( k , J ) 

35  continue 


(% 


* 


A 


f* 


cl  40  s  =  wa ( i ) *c dabs ( sum) 

40  s  =  wa  (  i  )  *ABS  ( sum) 

if  (f-  .  s  e .  s )  so  to  45 
f •  =  a 
imax  =  i 
45  continue 

c  test  for  a  lsor  ithmic  siri3ularitY 

if  <rn+F-  .  es.  rn)  so  to  105 
if  (J  .es.  imax)  so  to  60 

c  i nter chanse  rows  J  and  imax 

s  s  n  =  - 1 .  #  s  s  n 
cm 

do  50  k  =  1 7  n 

temp-  =  a(imax»k) 
a(imaxik)  =  a  ( J  »  k  ) 
a(Jik)  =  terriF 
50  continue 

wa ( imax )  =  wa ( J ) 

60  wa ( J )  =  i ma  x 

J  f  1  =  J  + 1 

if  (Jf-1  .st.  n)  so  to  70 

c  divide  by  F-ivot  element  u(J*J) 

temp  =  a  <  J  >  J  ) 
do  65  i  =  JF-l.n 

a  (  i  *  J  )  =  a  <  i  ,  J  )  / 1  erne 
65  continue 

70  continue 

cm  calculate  determinant 

do  72  i=l.n 

det  =  det*a(i»i) 

72  continue 
cm 

det  =  ssn*det 
cm 

75  if  (iJob  .  es.  1)  so  to  9005 
do  103  k  =  1 . m 

c  solve  ox  =  v  for  x 

iw  =  0 

do  90  i  =  l,n 
imax  =  wa ( i ) 
sum  =  b (imaxi k) 
b ( i ma  x  >  k  )  *  b ( i »  k  ) 
if  (iw  .es.  0)  so  to  B5 
iml  =  i-1 
do  80  J  =  iwiiml 

sum  =  sum-a ( i » J ) *b ( J » k ) 

80  continue 

so  to  88 

85  if  (t(l)  .ne.  zero  .or.  t(2)  .ne.  zero)  iw  =  i 

88  b(i»k)  =  sum 

90  continue 

c  solve  Iv  =  b  for  v 

n 1  *  n+1 
do  100  iw  =  1 » n 
i  =  n  1  -  i  w 
J  p  1  =  i  + 1 
sum  -  b  (  i  >  k ) 

if  ( J p 1  . s t •  n)  so  to  98 
do  95  J  =  J  f 1 .  n 

sum  *  s um-a ( i »  J ) # b ( J »  k ) 

95  continue 

98  b(i.k)  =  sum/a (i.i) 

100  continue 


103  continue 
so  to  9005 

105  ier  =  129 
9000  continue 

cal  l  uertst ( ier , 6hlestlc ) 
9005  return 
end 

subroutine  uertst  (ier, name) 
:  IMPLICIT  REAL*8(A-H,0-Z> 


algorithmic  singularity 


print  error 


c-uertst- 


■1  ibrarv 


function 

usase 

PARAMETERS 


lansuase 


latest  revision 


-  error  message  generation 

-  call  uertst ( ier » name) 

-  error  PARAMETER,  type  +  n  where 

tvpe=  128  implies  terminal  error 

64  implies  warning  with  fix 
32  implies  warning 

n  =  error  code  relevant  to  calling  routine 

-  input  vector  containing  the  name  of  the 

calling  routine  as  a  six  character  literal 
string. 

-  fortran 


-  January  18,  1974 


subroutine  uertst (ier, name) 
IMPLICIT  REAL*8(A-H,0-Z) 


dimension 

integer 

integer 

equivalence 


data 


ityp 


i  b  i  t 
printr 


data  printr  /  3/ 

i er2=i er 

if  (ier2  .se.  warn)  so  to  5 


ityp (5,4) , ibit (4) 
name (3) 

warn, war f , term, printr 

<ibit(l),warn),  (ibit(2),warf), ( i b i t ( 3 ) , t  er r ) 
/'warn' , ' ins  ' , ’ 

' warn ins (', 'with ',  '  fix',')  ’, 

' term' , ' inal  ’ ,  '  '  ,  '  '  ,  '  '  , 

' non-’ , ' def i ' , ' ned  ','  ','  '/, 

/  32,64, 128,0/ 

/  3/ 


ier 1=4 
so  to  20 
if  (ier2 

i er 1=3 
so  to  20 
if  (ier-2 

i er 1=2 
go  to  20 

ier 1=1 


It.  term)  90  to  10 


non-def ined 


termina 1 


It.  warf)  so  to  15 


warning (with  fix) 


warning 


extract 


i er2=i er2-i bi t ( i erl ) 

print  error  message 

write  (erintr, 25)  ( ityp ( i , ier 1 > , i=l , 5) , name, ier2, i er 

format ('  ***  i  m  5  1 {uertst)  ###  ',5a4,4x,a4,a2,4x,i2, 

*  '  ( ier  =  ’ , i3, ' ) ' ) 

return 


end 

subroutine  b  s  f  1  c  f 
C  IMPLICIT  REAL*8(A-H,0-Z> 

c####*  written  bv  J.  c.  wilev  univ.  of  texas  at  austin  Jan  3.976 

c - routine  computes  the  coefficients  of  the  b-selines  which  form 

c.  a  basis  over  the  set  of  knots*  r ,  with  repeated  knots 

c  at  the  end  points.  s(J,l,i)  j-power,  l-sesment,  i-seline. 


PARAMETER (NP  =21, 
l  neip=  1, 

3  n  s  i  p  =  1 , 

1  n  d  i  p  =  1 , 

5  ntip=  4, 

j  nsdp=101, 

7  nm2p=np-2,  nmlF=np-l ,  npl  F  =np  +  1 ,  np2p=np+2» 

7  nett5<ip=ne'ip*rie^f-,  nep=4*neqp+l , 

3  maxp=neqp*nsdF ,  nsap=16#np2p,  na  1  p=7tnetis sp , 

)  ndm2p=neqp#np2p ,  ndm3p=neqsqp*np2p , 

5  ndm23p=3#neqp ,  ndm33p=3#neqsqp) 


c  ommo  n  /  n c  1  /  nm2 ,  nm  1 ,  n ,  n  p  1 ,  n f-2 
common/r  c  1/r  (np  )  ,  rn  (np ) 
c  ommo n/ sac 1/sa (nsap) 

common/phiOl/F-Ol  (of)  /phi02/p02  (np)  /phi03/p03  (np) 
c  ommon/phi  1 1  /p  1 1  ( n f-  )  /  phi  12/p  12  (np )  /ph i  13/p  13  ( np ) 
common/phi21 /f-21  (nr)  /phi22/p22(np)  /phi23/p23  (np) 
c  ommon/phi i 1 /pi  1 (np) /phi i2/pi2 (np) 
common /phi  i3/pi3  (np)  /  phi  i4/pi4  (np) 
common/F  hi4/e (np2p) /phi5/f (np2p) 
common/bndvls/bcOp,  be  Of-  1 ,  be  1  f.-  1 ,  be  1  p ,  be  1 ,  bcO 
common/errc 1/err (np2p) 
common/serrc 1 /amxer 1 , amxer2, ermax 
c ommon/duml /d 1 (np2p) /dum2/d2 ( np2p ) 


commons  for  spline  intesrals 
c 

c####  ss  i  ( np2,  i  )  ,  dsi  ( riF-2,  7,  J  )  ,  ts  i  ( np] 
c####  i=n umber  of  sinsle  intesrals 
c####  J  =n umber  of  double  intesrals 
c####  k=number  of  triple  intesrals 


!,  49,  k ) 


common/ ss i c 1 /ss i (np2p, nsip) 
common/ dsic 1/dsi (op2p, 7, ndip) 
common/tsic 1/tsi (np2p, 49, ntip) 

c  ommo n / n s s i p / .  . i , n s J ( n s i p ) , ns v ( 2 , n s i p ) 
common /nds ip/ndi , ndJ (ndip), ndv(4, ndip) 
common/ritsip/nti  ,  ntJ  (ntip)  ,  ntv  (6,  ntip) 

equivalence  (sa,sa3) 
dimension  sa3 (4, 4, np2p ) 


c4=4. 00/ ( (r (2)-r (1) ) **4) 
i dx=16 

sa ( i dx )  =  -c  4 

sa(idx-l)  *  3.0*c4*r(2) 

sa  (  i dx-2)  =  -3.  0*c4*r  (2)  *r  (2) 


c4*r (2) *r (2) *r (2) 


sa ( i dx-3)  = 
i=2 

c4=4. 0/ ( (r  <3)-r (2) ) * (r (3)-r  ( 1) ) *#3> 
c3=4. 0/ ( (r(2)-r(3) )*(r(2)-r(l) )#*3> 
i dx=32 

sa ( i dx )  =  -c  4 

sa  <  idx-1 )  =  3. 0*c4#r  (3) 

sa  (  i dx-2)  =  -3. 0*c4*r (3) *r (3) 

sa  (  i  dx-3)  =  c4*r  (3)  *r  (3)  *r- (3) 

sa(idx-4)  =sa(idx)  -c3 

sa(idx-5)  =sa(idx-l)  +3. 0#c3#r (2) 
sa(idx-6)  =sa  <  i  dx-2)  -3.  0#c3#r- (2>  #r  <2> 
sa(idx-7)  =sa  <  idx-3)  +c3*r  (2)  #r- (2)  *r- (2) 

i=3 

c4=4.  0/  (  (r  (4)  -r  (2)  )  *  (r  (4)-r-  (3)  )  * ( r (4) -r ( 1 ) ) * *2) 
c3=4.0/((r(3)-r(2))*(r(3)-r(4))*(r(3)-r(l)  )**2> 
c 2=4.0/ < (r  <2)-r (3) >  * (r (2>-r (4) ) * (r <2)-r (1) ) **2) 
i  dx=48 

sa(idx)  =  -c4 

sa ( idx-1 )  =  +3.0*c4*r (4) 

sa  (  i  dx-2)  =  -3.  0*c4*r-  <4)  *p  (4) 

sa(idx-3)  =  c4*r (4) *r (4) *r (4) 

sa(idx-4)  =sa(idx)  -c3 

sa(idx-5)  =sa(idx-l)  +3. 0#c3#r (3) 
sa(idx-6)  =sa < i dx-2)  -3. 0*c3*r < 3) #r (3) 
sa(idx-7)  =sa (idx-3)  +c3#r (3) *r (3) #r <3) 

sa(idx-8)  =sa(idx-4)  -c2 

sa  (  idx-9)  =sa  (idx-5)  +3. 0*c2*r- (2) 
sa  (  i dx-10) =sa ( i dx-6)  -3. 0*c2#r (2) #r- (2) 
sa(idx-ll)=sa  ( idx-7)  +=2*r (2) *r (2) #r  (2)  • 

i=4, n-1 

do  10  i=4,nml 

ml=i-3 

m2=i-2 

m3=i-l 

m4=i 

m5=i+l 

c4=4. 0/  (  <r  (m5)-r  (ml) )  *  (r  (m5)-r  (m2)  )  *  <r  <m5) -r  (m3)  )  *  (r  (m5)-r  <rr.4>  )  ) 
c3=4. 0/  ( <r  (m4)-r-  (ml)  )  *  (r  (m4)-r  (m2) )  *  (r  (m4) -r  (m3)  )  *(r  (m4)-.r  (m5) )  ) 
c2=4. 0/  (  < r  ( m3 )  - r  (ml )  )  *  (r  (m3)  -r  (m2)  )  *  (r  (m3)  -r  (m4)  )  *  (r  (m3)-r  (m5)  )  ) 
c  1=4.0/  (  ( r  (m2) -r  (ml)  )  *  (r  <m2)-r-  (m3)  )  *(r  (m2)-r-  (m4)  )  *  (r  (m2)-r  (m5>  )  ) 
idx=16*i 

sa(idx)  =  -c4 

sa (idx-1)  =  +3. 0#c4#r (m5) 

sa(idx-2)  =  -3.  0#c4*r  (m5>  $r-  (m5) 


sa(idx-3)  = 


c4*r  (m5)  #r  (m5)  *r  (m5) 


sa(idx-4)  =sa(idx)  -c3 

sa(idx-S)  =sa(idx-l)  +3. 0*c3*r (m4) 

sa(idx-6)  =sa  (  i  dx-2)  -3. 0#c3#r (m4) #r (m4) 

sa(idx-7)  =sa(idx-3)  +c3*r  (m4)  #r  (m4)  #r  <rn4) 

sa(idx-8)  =sa(idx-4)  -c  2 

sa(idx-9)  =sa(idx-5)  +3. 0*c2*r (m3) 

sa ( i dx-10) =sa ( i dx-6)  -3. 0#c2#r (m3) #r (m3) 

sa  ( i  dx-1 1 )  =  sa  (  i  dx-7)  +c2#r (m3) #r (m3) *r (m3) 

sa ( idx-12) =  sa ( i dx-B)  -cl 

sa(idx-13)=sa(idx-9)  +3. 0*c 1 #r (m2) 

sa(idx-14)=sa(idx-10)-3.0#cl#r(m2)H<r(m2) 

sa  ( idx-15)  =sa  ( idx-1 1 )  +c  1  *r  (m2)  *r  (m2)  *r  (m2) 

10  continue 
-i=n 

ml=n-3 
m2=n— 2 
m3=n-l 
m4=n 

c4=12. 0/ ( ( r ( m4 ) - r (ml) ) # (r (m4)-r (m2) ) *  (r (m4)-r (m3) ) ) 
c3=-4. 0/ ( (r (m4)-r (m2) ) * (r (m4)-r (m3) ) * (r (m4)-r (ml) ) **2) 

2  -4.0/  (  (r  <m4)  -r  (ml) )  *  (r-  (m4)-r  (m3) )  *  (r  (m4)-r  (m2)  )  **2) 

3  -4.0/ ( (r (m4)-r (ml) ) * (r (m4)-r (m2) )  *  (r (m4)-r (m3) ) **2) 
c2=  4.0/ ( (r (m3) -r (ml) ) *  (r (m3)-r (m2) >  * ( r (m3)-r (m4) ) #*2) 
cl=  4.0/  (  (r  (m2)  -r  (ml)  )  *  <r  (m2)  -r  (m3)  )  *.(r  (m2)-r  (m4>  )  **2) 
idx=16#n 

sa(idx-4)  =  -c  3 

sa(idx-5)  =c4  +3. 0*c3*r (m4) 

sa(idx-6)  =  -2.0#c4  -3. 0*c3#r (m4) (m4) 

sa  (  i  dx-7)  =c4*r  (m4)  *r  (m4)  +c3#r-  <m4)  #r  (m4)  #r  (m4) 

sa(idx-8)  =sa(idx-4)  -c2 

sa(idx-9)  =sa(idx-5)  +3. 0*c2#r (m3) 

sa  < i dx-10) =sa ( i dx-6)  -3. 0*c2#r (m3) #r (m3) 

sa (idx-1 1 ) =sa ( idx-7)  +c2*r (m3) *r (m3 ) *r (m3) 

sa ( i dx-12) ^sa < i dx-8)  -cl 

sa(idx-13)=sa(idx-9)  +3. 0*c 1 *r (m2) 

sa(idx-14)=sa<idx-10)-3.0*cl#r(m2)#r(m2) 

sa  (idx-15)  =sa(idx-ll)  +c  1  *r-  (m2)  #r  (m2)  #r  (m2) 

-i=n+l 

c4=  12.0/ (  (r  (m4)-r  (iti2)  )  *  (r  <m4) -r  (m3) ) ) 

c3=-12. 0/ ( (r (m4)-r (m3) )  * (r (m4) -r (m2) ) **2> 

2  -12.0/  (  (r  (rr.4)-r  (m2) )  *  (r-  (m4) -r  (m3) )  **2> 

c2=  4. 0/ < ( r (m4>  — r (m3) >  * (r (m4) -r (m2) ) **3) 

2  +4.0/  ( (r  (rr.4)-r  (m2) )  **2  *  (r  (m4> -r- (m3)  >  **2> 

3  +4.  0/ ( (r (m4) — r (m2) )  * ( r (m4) -r (m3) ) **3) 

c 1=  4. 0/ ( (r (m3) -r (m2) )  # ( r (m3) -r <m4) ) ##3) 

idx=16*(n+l) 

sa(idx-8)  *  -c2 

sa(idx-9)  *  c3  +3. 0*c2*r (m4) 

sa  ( i  dx-10)  =-c4  -2. 0*c3#r  (m4)  -3. 0*c2*r  (m4)  * r-  (m4) 

sa ( idx-1 1 ) =  c4#r  (m4)  +c3#r  <m4)  #r  (m4)  +c2*r  (m4)  #r  (m4)  *r  (m4) 

sa ( idx-12) =sa ( i dx-8)  -cl 

sa(idv-13)=sa(idx-9)  +3. 0#c 1 #r (m3) 

sa ( idx-14) =sa ( i dx-10) -3. 0*c l*r (m3) *r (m3) 

sa ( idx-15) =sa ( i dx-1 1 )  +c l*r (m3) *r (m3) *r (m3) 

-i=n+2 

c4=  +4.0/ (r- (m4) -r  (m3) ) 
c3=-12. 0/ (r (m4> -r (m3) ) **2 
c2=+12.0/(r(m4)-r (m3) )**3 
cl=  -4.0/  (r-  (m4)  -r  (m3)  )  **4 
idx=16*  <r.+2) 

sa(idx-12)*  -cl 

sa(idx-13)=  c2  +3. 0*c 1 #r (m4) 

sa(idx-14)=  -c3  -2.  0#c2!t:r  (m4)  -3. 0*c  1  #r- (m4)  #r  (m4) 


sa  (idx-15)  =c4+-c3#r  (m4)  +c2#r  (m4 )  *r  (m4)  +c  1  *r  (m4 )  *r  ( m4)  #r  (m4 ) 

return 

end 

subroutine  depse(fn,c) 

C  IMPLICIT  REAL*8(A-H,a-Z) 

c  #  #  #  #  *  written  bY  J.  c.  wileY  univ.  of  texas  at  austin  Jan  1976 
dimension  fn(5),c(l) 


c 

c#### 

c#### 

c#### 

c#### 


PARAMETER (NP  =21, 

2  ne«*p=  1, 

3  nsiF=  1, 

1  ridip=  1, 

5  ntip=  4, 

3  ri3dF=101, 

y  nm2p=np-2, nml p=rip-l , np 1 p=np+l , riF  2p=np+2, 

y  neasqp=nesplcnesp ,  nep=4*neap+l , 

3  maxp=neap#nsdp, nsap=16*np2p, nalp=7#neasap, 

1  ndm2p=neap*np2p,  ndm3p=neqsqp#np2p , 

5  ndm23p=3*neap,  ndm33p=3Htneasap> 


common/nc  l/nm2»  nml ,  n»  nF  l ,  riF-2 
c  ommo  n / r  c 1 / r ( n  p ) , r  n ( n  p ) 
common /sac 1/sa (nsar) 

common/phiOl  /  f-01  ( np) /Phi02/p02 (np )  /Phi03/p03  <  np ) 
common /phi  Il/pll(np)/phil2/pl2(rip)/phil3/pl3(np) 
common/phi21/p21 (np) /p h i 22/f22 ( np ) /phi23/p23 (np) 
common/ Phi i  1  /  pi  1 (np)/phii2/pi2(np) 
common/ phi i3/pi3(np) /phi i4/Pi4 (np) 
common/phi4/e(np2p) /Phi5/f (np2p) 
common/bridvls/bcOp,  bcOp  1 ,  be  Ip  1 ,  be  1  p,  be  1 ,  bcO 
common/errc 1/err (np2p) 
common/serre  1/amxerl,  amxer-2,  ermax 
common/duml/dl (np2p) /dum2/d2(np2p) 

commons  for  spline  integrals 

ssi (np2» i) ,dsi (np2,7,J) ,tsi (np2,49, k ) 
i=number  of  sinsle  integrals 
J  =number  of  double  intesrals 
k=number  of  triple  intesrals 

common/ssic l/ssi (np2p,  ns  ip) 
common/dsicl/dsi (np2p,7,  ndip) 
common/tsic 1 /tsi ( np2p,  49,  nti p) 

c  ommo n / n s s i p / n s i , n s  J <  ns i p ) , n s  v  < 2 , n s i p ) 
common/ndsiF  /ndi , ndJ (ndip) , ndv  <4, ndip) 
common/ntsip/nti ,  ntJ  (ntip) ,  ntv  (6,  ntip) 

equivalence  (sa,sa3) 
dimension  sa3 (4, 4, op2p) 


nm3=n-3 

if (abs(fn(l) > . st.abs <10. *fn (2) > > 
■find  derivative  of  fen  at  end  Pts. 
dll»<fn<2)-fn(l) ) / (r (2)-r <1) ) 
dl2"(fn(3)-fn(2) )/(r(3)-r(2) ) 
dl3=(fn(4)-fn(3))/(r(4)-r(3) ) 
d21=(dl2-dl 1 ) / (r (3>-r (1) ) 


d22=  <d 13-d 12) /  <r  (4) -r (2) ) 
d31=  ( d22-cJ21 )  /  <r  (4)  -r  (1  >  ) 
f  p0=dll-d21*r- <2)+d31*r  (2)  *r  (3) 


* 


* 


I  ». 

»  • 


I 


* 


I 

I 


dl  1  =  ( f  n  ( nm2 ) -f  n  <  nm3>  )  /  (r  (nm2)-r  (  nrn3)  ) 
d  1 2=  ( f  n  ( nm  1 )  —  f  r.  ( nn.2)  )  /  ( r  <  nm  1 )  -r  <  nrr.2)  ) 
d 1 3=  <  f  n ( n ) -f  n  <  nml ) ) / ( r ( n  >  -r  ( nm 1 ) > 
d21  = (dl2-dl 1 ) / (r (nml > -r < nm3) ) 
d22=  < d  13-d  12)  /  (r  (r.)-r  <nm2>  ) 
d31  =  (d22-d21 )  /  <r  (n)  -r  (r.m3)  ) 

f>l=+r  <  n >  *  (dl3+(r  (n)-r  (nml )  )  *  <d22+d31*  ( r-  (n)-r  (nm2)  )  )  ) 

c - compute  e  and  f. 

c ( 1 ) =  f  n ( 1) /fOI (1 ) 

22  continue 

f  (2)  =  <f  pO-p  1 1  (1)*c(1))/p12(1) 
e (2) =0. 0 

e  (3)  =-f-03  (2)  /f02  <2) 
f  (3)  =  (fn  (2)  -p-01  (2)  *f  (2)  )  /f-02<2) 
do  10  i =4 , n 

e  <  i  )  =1 . 0/  <  f-01  (i-l)  *e(i-l)+F  02(i-l)  ) 
f  (  i  )  =  (f  n  (  i-l )  -f-01  (i-l)*f(i-l))*e(i) 
e(i)=-p03(i-l) *e<i) 

10  continue 

c - compute  c. 

c  ( np2)  =f  n  ( n)  /f-03  (n) 
c(npl)  =  (fpl-pl3(n)*c  (r.F-2)  )  /pl2(r.) 
do  12  i=2\  n 
j =np2-i 

c  ( J ) =e ( J ) # c ( J  + 1 ) +f ( j ) 

12  continue 
return 
20  continue 

dll=(fn(3)-fn(2> ) / ( r <3) -r (2) ) 
dl2=(fn (4) — f n (3) ) / (r <4)-r (3) ) 
dl3=(fn  <5)~fr.  (4)  )  /  (r  <5>-r  (4)  > 
d21= (dl2-dl 1 ) / (r (4) -r (2) ) 
d22= (dl3-dl2) / (r (5)-r (3) ) 
d31=(d22-d21) / (r  <5) -r (2) ) 
f nO=fn (2)  +  (r  <l)-r (2) ) # (dl l+(r (l)-r (3) ) *(d21 
2  +<r(l)-r(4) )*  d31)) 

fpO=dll  +  d2U  (  (r  ( l)-r  (3)  )  +  (r  <  1 ) -r  (2)  )  >+d31*  (  ( r  ( 1 )  -r  (3)  )  * 

2  <r  < 1) -r (4) )  +  (r (1 ) -r (2) ) # (r (l)-r  (4)  )  + 

3  <r  <  1 )  -r  (2)  )  #  (r  ( 1 )  -r  (3) ) ) 
c ( 1 ) =f nO/FOl (1) 

so  to  22 
end 

subroutine  deset (Jh) 

C  IMPLICIT  REAL*8<A-H,0-Z) 

ctttt*  written  bv  J.  c.  wilev  univ.  of  texas  at  austin  Jan 

c 

c 

PARAMETER (NP  =21, 

2  nesp=  1, 

3  nsiF=  1, 

4  ndie=  1, 

5  ntip=  4, 

6  ri9dp=101, 

7  nm2p=np-2,  nml  p=np-l ,  nr  1  p=hp+1  ,  n F-2p=nr+2, 

7  nesssp=nesF  *nesp , neF  =4*nesp  +  l , 

3  maxp=nesf-*n9dp,  nsap=16*rip2p,  nalp=7#nesssF  , 

4  ndm2p=nesp*np2p »  ndm3p=ness«ip*np2p , 

5  ndm23F  =3*ne*iF  ,  ndm33r ^SfcnesssF  ) 


1976 


c  ommon/  ric  1  /  nm2 ,  nml  >  n>  hf  1 ,  n f  2 
c  ommon/ rc  1  /  r  (  n f  )  ,  rn  (of) 
c  ofTimon/ sac  1  / ss  (  n s a f  ) 

c  ommon/ PhiOl /pOI ( np ) /ph 102/ p02  <  np )  /f  h i03/F  03  <  hf  ) 
c  ommo n  / p h  i  1 1  / f  1 1  ( n f  )  / p  h i 1 2/ f 1 2 ( n  r )  /  f  •  h i 1 3/ f 1 3 ( op ) 
common/ F-hi21/F-21  <r.F  >  /ph i22/F-22  <  of-  )  /f  h  123/f  23  < hf  ) 
c  ommo n/ F-h i  i  1  / f-  i  1  (  of  )  / f  h  i  i 2/f  1 2  (  of  ) 
c  ommo n/ F-h i i3/pi3  <nF ) /phi i4/p j  4 (of ) 
common/phi4/e (np2p) /phi5/f (np2p) 
common/bndvls/bcOp,  bcOr  1 ,  t>c  1f-1  ,  be  lF  >bcl>  bcO 
c  ommon/erre 1 /err ( np2p ) 
c ommon/serre  1  /amxer  1  .  amxer 2i  er  max 
c  ommo n /d urn  1  /  d  1  ( np2p )  /dum2/d2  ( np2p ) 
c 

c  commons  for  spline  integrals 

c 

c####  ss i ( of  2i i ) »  ds i ( nr2>  7, j ) >  ts i ( np25  49 >  k ) 

c####  i=number  of  single  integrals 

c####  J=number-  of  double  integrals 

c####  k-number  of  triple  integrals 

c 

c ommon/ss i c 1 /ss i < np2e  >  ns i p ) 
common/ dsi c 1/d si (np2p,  7 >  ndiF ) 
c  ommo  n  / 1  s  i  c  1  / 1  s  i  ( n p 2f  i  49  f  n t  i  F- ) 
c 

c  ommo  n / n  s  s i p / n s i * n  s J ( ns i p ) »nsv(2»nsip) 
c ommon/nds i F/nd  i  .  ndJ  (ndiF  )  .  ndv  (4.  ndiF  ) 
common/ritsip/riti  .  ntJ  (ntip) ,  ntv  (6.  ntip) 
c 

equivalence  <sa,sa3) 
dimension  sa3 (4, 4, np2p) 
c 
c 

data  p  i  1  ( 1 )  i  p  i2  ( 1 )  »  pi.3  ( 1 )  i  F-i4  ( 1 )  /4#0.  0/ 

c - subroutine  sets  up  values  of  splines  at  the  knots. 

do  10  i=l,n 

call  spvI  ( i  i  r  (  i  ) »  f-01  ( i ) » 1 ) 
call  spvI  <  i  +  1 .  r  ( i )  *  f-02  ( i ) »  1 ) 
call  s F' v  1  ( i  +2 1  r  ( i  )  i  f-03 ( i ) » 1 ) 
call  spvlp(i»r(i).pll(i).l) 
call  SPvlp<i  +  l7r<i)»F-12(i)>l) 
call  s p v 1 f  < i +2 *  r ( i ) »  p 1 3 ( i ) » 1 ) 
call  sp v  1  F-p  ( i  »  r  (  i  )  >  f-21  (  i  )  7  1 ) 
call  sp vl pp ( i  +  1 *  r ( i ) »  p22 ( i )  ,  1) 
call  spvIpp  <i+2, r (i )  ,  p23(i )  ,  1) 

10  continue 
p02 ( 1 ) =0. 0 
f-03  <  1 )  =0.  0 
pOI <  n ) =0 . 0 
f-02  (n)  =0.0 
f-13  ( 1 )  =0. 0 
p 1 1 ( n ) =0 . 0 
do  12  i=2.n 
F-i  1  (  i  )  =gaus6  <  i-1 .  i  ,  j h ) 

Pi2(i)=saus6(i7  i  J  h) 
p13 ( i ) =saus6 ( i  +  1 . i . J  h) 

Pi4 ( i ) =saus6 ( i+2» i » J  h) 

12  continue 

be  Of-  =  f  11  (1) 
t-cOrl  =  f  1 2  ( 1 ) 
belFl  =  p 1 2 ( n ) 


r  1 3  ( n  ) 
F-  03  ( n ) 
f01  (1 ) 


ti  C  1  F  = 

be  1 

beQ  = 
retu rn 


end 

function  ea  us  10  (  f  1  ,  p2»  p3,  n  ,  a  ,  b ) 

C  IMPLICIT  REAL*8 (A-H, 0-Z) 

c***#*  written  bv  J.  c.  wilev  univ.  of  texas  at  austin  Jan  1976 
dimension  p  1  (5)  >f-2(5>  >  p3  (5) 

datawl,w2,w3/. 467913934572691, .360761573048139, . 171324492379170/ 
datax 1 , x2, x3/. 238619186083197, . 661209386466265, . 932469514203152/ 
fen  (x)  =  (f-1  (i)+x*  (f  1  (2)+x¥  (el  <3>+x#  (el  <4)+x*rl  (5)  )  >  )  )  * 

2  <f-2(1  )+x*  <p2(2)  +x*  (f-2(3)+x*  (f-2  (4) +x*p2  (5)  )  )  )  )  * 

3  <p3(l)+x*(r3(2>+x*(p3(3>+x#(p3(4)+x*e3(5) ) ) ) )* 

4  (x*#n) 
rd=0.5* (b-a) 
rp=Q. 5* (b+a) 

sauslO  =r  d¥  ( w3#  ( f  cn ( rp— r  d¥x3) +f  c  n ( rp+rd*x3) ) 

2  +w2# (fen (rp-r d*x2) +f cn (rr+r d*x2) ) 

3  +wl # (f c  n ( rp-r  d#x 1 ) +f c  n ( r F  +r  d*x 1 ) ) > 


return 

end 

function  saus6 (k > i , J h ) 

C  IMPLICIT  REAL*8<A-H,0-Z> 

c#**#¥  written  bv  J.  c.  wilev  univ.  of  texas  at  austin  Jan  1976 


c 

c 


PARAMETER (NP  =21, 

2  nesp=  1, 

3  n  s  i  r=  1 , 

4  ndip=  1, 

5  ntip=  4, 

6  n  s  d  p  =  1 0 1 » 

7  nm2p=np-2,  nm).p=np-l ,  riFlp=nP+l ,  np2p=np+2, 

7  neis«ip=nesp*rie^p.  nep=4i|irie«iP+l , 

3  maxp=nee)p#nsdp,  nsap=16*np2p»  nalp=7*negs*p, 

4  ndm2p=nesp*np2p»  ndm3F=negssp #np2p  , 

5  ndm23p=3#nesF , ndm33p=3#nesssp) 


e 

c 

c 

common /nc l/nm2, nml , n , npl , np2 
common/rc 1/r (np) , rn (np) 
c  ommcin/sac  1 /sa  ( nsap ) 

c  ommon/ph iOl /pOI <  n  p ) /ph i02/p02 ( np ) /phi03/p03 ( np ) 
c  ommon/ph ill/pll(np)/phil2/pl2(np)/phil3/pl3(np) 
common/ phi 21 /p21 (np) /phi22/p22 (np) /phi23/p23 (np) 
c  ommo n  / F- h  i  i  1  /  p  i  1  ( n p  )  / F- h  i  i  2/  p  i  2  ( n  p  ) 
c  ommo  n / p  h i i 3 / p i 3 ( n  f ) / p  h i i 4 / p i 4 ( n  p ) 
c  ommon/ph i 4/e ( np2p ) /Phi5/f ( np2p ) 
common/bndvls/bcOp, bcOpl , be IpI , be  Ip, be  1 , bcO 
common/errc 1/err (np2r) 
common/serrc  1  /amxer  1 ,  amxer-2,  ermax 
c ommon/duml / dl  ( nr2p )  /durn2/d2  ( np2p ) 
c 

c  commons  for  spline  intesrals 

c 

c####  ssi  (  nF  2,  i  )  ,  ds  i  (riF 2,  7,  J  )  ,  tsi  ( np2,  49, 1: ) 
c####  i=number  of  single  integrals 
c.####  J=riumber  of  double  intesrals 
c####  k=n umber  of  triple  integrals 
c 

c  ommon/ss i e  1  /  ssi  inr  2p ,  ns,  i f  ) 


c  orrirrion  /  6  s  i c 1 / ds  i  ( np2r , 7 , nd i f ) 
common /t-sicl/tsi  ( n f 2p ,  49,  rt t  i  f- ) 
c 

c  ommort/ns s i p/ns i  ,  nsj  <  ns  i  p  )  ,  ns v  (2,  ns i  p ) 
c  nmmon/nds i  f  /  nd  i  ,  ndJ  ( rid  i  p  )  ,  nd  v  (4,  nd  i  p  ) 
c  ommon /ntsip/nti, n t J ( n t i p ) ,  n  t  v  ( 6 ,  n t i p > 
c  . 

equivalence  (sa,sa3) 
dimension  sa3 (4, 4, np2p) 


c 

c 

da taw 1,  w2,w3/.  46791 3934572691, .360761573048139, . 171324492379170/ 
dataxl , x2, x3/. 238619186083197, . 661209386466265, . 932469514203152/ 
dimension  y(6),x(6) 
r  d=0. 5# ( r ( i ) -r < i-1 ) ) 
r f=0.  5#  ( r  <  i  ) +r  (  i-1 )  ) 
x ( 1 ) =rp-rd#x3 

x  (2>=r-p-r  d*x2 
x  (3)  =r-F  -rd#xl 
x (4) =r  p+r  d*xl 
x  (5)  =rp+r d#x2 
x (6) =r F+r d#x3 
call  spv1U:,x(1),y<1),6) 

saus6=rd*  (w3# ( x ( 1 ) **J h# y ( 1 ) +x  <6) #*J h# v  <6) ) 

2  +w2* (x (2) **ih*v (2) +x (5) **Jh#v (5) ) 

3  +wl  # (x (3) *#J  h#v (3) +x (4) #*J  h*v  (4)  )  ) 


return 

end 

subroutine  arid (er, isw) 

C  IMPLICIT  REAL*8(A-H,0-Z) 

•  c***#*  written  bY  J .  c.  wilev  univ.  of  texas  at  austin 
c 
c 


3 

4 

5 

6 
7 
7 

3 

4 


PARAMETER (NP  =21, 
neqp=  1 , 
nsip=  1, 
ndip*  1* 
ntip=  4, 
nadp=101 , 

nm2p=np-2,  nml.F  =np-l  ■»  rip  1  p=np+l ,  np2p=np+2. 


neqsqpsneqplneqp, nep=4#neap+l , 
fnaxp=ne*»p#nsdp ,  nsap=16#np2p,  na  1  p=7#ne«»ssp, 
ndm2p=neap)Knp2p ,  ndm3p=ne«i5<qp#np2p , 
ndm23p=3#ne<qp,  ndm33p=3#ne<qs,qp ) 


Jan  1 976 


c 

c 

c 

c  ommon /nc  l./nm2,  nml ,  n ,  rtP  1 ,  ne2 
m.  common/rcl/r (np)  ,  rn (np) 

c  ommon/sac l/sa (nsap) 

common/phiOl/F-Ol  (np )  /phi02/p02 (np)  /phi03/p03  ( np) 
c ommon /phi  1 1  /  p1  1  (r,p)  /phi  12/p  12  (np)  /  phi  13/p  13  (np) 
c  o  mm  o  n / p  h i 2 1 / p  2 1 ( n  p ) / p  h i 22 / p  2 2 ( n  p ) / f  h i 2 3 / p  2 3 ( n  p ) 
c  ommon /phi i 1 /p i 1 ( nr) /phi i2/p i2 ( np) 

•  common/ phi i3/pi3(np)/phii4/pi4(np) 

common/phi4/e(np2p) /phi5/f (ne2p) 
common/bndvl s/bcOp , be Op  1 , be  1 p 1 , be  Ip, be  1 , bcO 
common/errc 1 /err ( np2p ) 
common/serrc  1/amxer- 1 ,  amxer-2,  ermax 
c ommon/ dum 1 /dl ( nF  2p ) /dum2/d2 (np2p ) 


c 

c 


commons  for  spline  intesrals 


c  ####  s  s  i  ( n f-2 ,  i  >  ,  d  s  i  (  n f  2 ,  7 ,  J  >  ,  t  s  i  ( n f  2 ,  49 ,  k  > 

c####  i=number  of  sinsle  intesrsls 

c####  J  =number  of  double  integral? 

c####  k=number  of  trifle  integrals 

c. 

c  otnmori/  ssicl/ssi  (  nr  2f-  7  nsi  f  ) 
c  o  mm  o  n  /  d  s  i  c  1  /  d  s  i  ( n  p  2  f  7  7  7  n  d  i  f  > 
c  o mm 0 n  / 1 s  i  c  1  / t  s  i  ( hp2p , 49?  nt i p ) 
c 

c  oiTifiio n /  nssi  f  / ns i  7  n  s  J  (nsir)  , ns v (2, ns i f ) 
c  ornmo n/nds i p/nd  i  7  ndJ  ( nd  i  p  )  7  rid  v  ( 4 7  nd  i p  > 
c  ommon/nts i p/nt i , nt J  ( nt  i  p )  7  n t  v  <67  nt  i p ) 
c 

equivalence  (sa7ss3) 
dimension  sa3 (4, 4? np2p ) 
c 
c 

dimension  er(l) 

c - 9rid  sets  up  srid  spacing. 

c  isw=l7  uniform  srid. 

c  isw=2, sfsc ins  based  on  err  fen. 

c  1 sw=3,  >  >>>>>>> 

c  i sw=4i un i f orm  except  end  ets. 

so  to  (301,802,803,804) , isw 
c - sets  up  uniform  mesh. 

801  do  10  i=l,n 

10  er ( i )  =  ( i-1 . 0) / (n-1 . 0) 
return 

c - this  section  chose?  a  new  set  of  knots  based  on  th 

c  function  er. 
c - 

c - note:  er  on  exit  contains  new  x. 

802  continue 
emax=0. O 

do  19  i~2 , nml 

19  emax-MAX (emax, er ( i ) ) 
emin=. 001  * emax 

do  20  i=2, nml 

er ( i ) =MAX ( emi n, er ( i ) ) 

er  (  i  )  =er-  (  i  )  *#0.  25 

20  continue 

er  ( 1 ) =er (2) 
er.(  n)  =er  ( nml ) 
sum-0. 0 
do  21  i=l,nml 

21  sum=sum+0. 5#  ( er  <  i )  +er-  (i+l))*(r(i+l)-r(i)) 
sum=sum/nml 


c - compute  partition. 

r  n ( 1 >  =0 . 0 
k  =  l 

tot=0. 0 
if  1 9*1 

do  22  i =2 7  nml 
soal= ( i-1 ) *sum 
25'  del  =  (r  < k  +  1  >  -r  ( k  >  ) 

add=del* (er <k+l >+er (k) ) *0.5 
if (add+tot. It. soal )  90  to  23 
rn ( i >  =r (k )  +  ( 90a  1-tot >  *deJ /add 
90  to  22 
23  tot*tot+add 
k*k  +  l 
90  to  25 


A 


A 


A 


22  continue 
rn ( 1 )=0. 0 
r  n  <  n )  =  1 . 0 
do  33  i=l,n 
e  r  (  i  )  =  r  n  ( i  > 

33  continue 
r etur n 

803  continue 

c - map  er  to  seal  ins  form. 

er  ( 1 ) =er (2) 
er ( n ) =er ( nml ) 
do  40  i = 1 i n 

if < er ( i ) . st . amxer2)  so  to  41 
if (er ( i ) . st . amxer 1 ) so  to  42 
i f < er ( i ) . st . 0. 2*amxer 1 )  so  to  43 
er ( i ) =1 . 25 
so  to  40 

41  er ( i ) =0. 60 
so  to  40 

42  er ( i ) =0. 75 
so  to  40 

43  er(i)=1.00 

40  continue 

c - compute  new  r. 

r i n=0. 0 
do  50  i si, nml 

ripln=r  in+  (r  (  i  +  1 )  —  r  <  i  )  >  #0. 5#  (er  (  i  )+er-  (  i  +  1 )  ) 
er  < i >  »r i n 
r i n=r i p 1 n 

50  continue 

er (n)«riPln 

c - rescale. 

do  51  i  ■  1 1  n 

51  er  (  i  )  =  er  (  i  ) /er-  (n) 

c - check  er  for  minimun  sracins. 

do  54  i=2>n 

if (er < i ) -er < i-1 > . st. .01)  so  to  54 
er  <  i )  =er-  ( i-1 ) +. 01 

54  continue 

do  55  i=2 . n 

55  er < i ) =er < i > /er (n) 
return 

804  continue 
nm4=n-4 
dlr=l . / (nm4-l ) 
er ( 1 ) =0. 0 

er (2) =. 5*dlr 
er (3) =d 1 r 
er (4) =1 . 5*dlr 
do  87  i=5.nm4 
er ( i ) = ( i-3) #d 1 r 

87  continue 

er ( n-3) = (nm4-2. 5) *dlr 
er ( n-2)  = (nm4-2) #d  lr 
er  (n-1 )  *  (nm4-l  .5)  ♦dir¬ 
er  ( n)  *1 . 
return 
end 

subroutine  p  o  p ( a »  p » 1 ,  rn  n  s ) 

C  IMPLICIT  REAL*8(A-H,0-Z> 

cf****  written  bv  J.  c.  wileY  uriiv.  of  texas  at  austin. 
dimension  p(5)»a(4) 
imax=4+n-l 


Jan  1976 


do  10  1=1,5 

10  p  ( i ) =0 . 0 

vi  m  i  n=ma  xO ( 1 + l ~n ,  1) 
ns=maxO  ( 1-jmi  n  +  1 , 0) 
i f ( j  m i n . s  t . 4 )  retur n 
do  11  J=Jmin»4 
f ac=l . 0 

i f ( 1 . es . 0)  so  to  12 
do  13  l 1=1 , 1 
13  f ac= ( J +n-l-l+i 1 ) #f ac 
12  f-  ( j  +ns-l )  =f ac #a  (J  ) 

11  continue 
return 
end 

subroutine  rear- dr  ( u ,  nk  ,  r  new) 
C  IMPLICIT  REALMS (A-H,D-Z) 

written  by  J.  c.  wilev  u 
c - reordr  interpolates  u  from  cur 


u  n  l  v . 
current 


of  texas  at 
s r i d , r ,  to  i 


a  u  s  t  i  ri  J  a  n 
i  ew  s  r  i  d . 


1976 


PARAMETER (NP  = 
l  n  e  s  f  = 

3  nsiF= 

\  n  d  i  f  = 

5  ntiP= 


nsdF  =101 , 

nm2p=np-2,  nml f  =nr-l ,  n f-  1  p=riF  +  1 ,  of  2f  =  np+2, 
nesssF  =nesF  #nesp , nep=4#negp+l , 
maxF-=nesF  CnsdF  ,  nsaF=16*np2p,  nal p=7#nensgF , 
ndrri2p=neap # n p2f-  ,  ndm3p=negs«»F  t np2e , 
ndm23p=3#negp,  ndm33p=3#  negssp ) 


c 

c#### 

c#### 

c#### 


dimension  rnew(l)  ,  u  <  n f-  »  ok ) 


c omrnon/nc  1  /nm2,  nml , n, npl , of  2 
common /r c 1 /r (np ) *  rn ( nr ) 
common /sac 1/sa (neap) 

c ommon/ph iOl /f 01  (nr)  /phi02/p02  (np)  /  phi03/r03 (nr) 
c  ommon/F  hi 1 l/pll(np)/phil2/pl2(np)/phil3/pl3(np> 
common/phi21/p21 (nr) /phi22/p22 (np) /phi23/p23 (nr) 
commori/F  hi  i  1  /Fi  1  (nr)/phii2/pi2(np) 
comirion/F  hi  i3/Pi3  (np)  /  phi i4/pi4  (nr) 
common/F-hi4/e  (hp2p)  /f  hi5/f  (nr2r) 
common/bndvls/bcOp,  bcOrl ,  be  1f-1  ,  be  Ip,  be  1 ,  bcO 
common/errc 1/err (hp2p) 
c  ommon/serre 1 /amxer 1 , amxer2, errnax 
common/duml/dl (np2r) /dum2/d2 (nr2p) 

commons  for  srline  integrals 

ssi (nr2, i  > , dsi (np2, 7, J ) , tsi (nr2, 49, k ) 
i=number  of  single  integrals 
J=number  of  double  integrals 
k=number  of  triple  integrals 

common/ssic 1/ssi (np2p» nsip) 
c  ommon/ds i c 1 /ds i (nr2p , 7*  ndi p) 
c  ommon/tsic 1 /tsi <  np2r , 49, nt ir ) 

common/nss  ir/ns i ,  ns.i  (nsi  f-  )  ,  ns  v  (2,  ns  ir ) 
common/ndsi  p/nd  i ,  rid  j  (ndip)  »ndv(4»ndir) 
c  ommon/nts l p/nt i , nt J (nt i r) , nt v (6, nt i p ) 


e^u  i  va  1  enc  e  (sa,sa3) 
dimension  sa3  (4, 4,  iip2p) 
dimension  u  i  ( 1 )  ,  sr  ( 1 ) 
c 
c 

do  12  k=l i nk 
call  depse(udik).ui) 
do  13  i  =  1 ,  n 
13  u ( i , k ) =0 . 0 
do  12  1=1, nr2 
call  s  p  v  1  ( 1 ,  r  n  ew ,  s  p »  n ) 
do  12  i=l,n 

u  (  i  , k ) =u ( i , k ) +ui ( 1  >  *sp  (  i ) 

12  continue 
return 
end 

subroutine  repse(c,v) 

C  IMPLICIT  REALMS ( A-H, 0-Z > 

c#####  written  fc<Y  J.  c.  wiley  univ.  of  texas  at  austin  Jan 
dimension  c(l)»v(l> 
c 

PARAMETER (NP  =21, 

2  n  enp=  1, 

3  nsip=  1, 

4  nd i p=  1 , 

5  ntip=  4, 

6  nsdp=101» 

7  nm2p=np-2, nmlF=np-l , nplp=np+l , nr2r=np+2, 

7  ne«tssp=ne‘:ip*ne>ip,  nep=4*ne«»p+l , 

3  maxp=nec»p*risdp ,  nsap=16# np2p ,  na  1  p=7#nee»ssp , 

4  ndm2p=neqF#np2p , ndm3p=nessgp#np2p> 

5  ndm23p=3Htne«»p ,  ndm33p=3#nec»sgp) 
c 
c 
c 
c 


c 

c 

c 

c#### 

c#### 

c#### 

c*«*« 

c 


c  ommo n  /  n c  1  /  nm2 ,  nm  1 ,  n ,  ri p  1 ,  n p 2 
common/r  c  1/r  (rip)  ,  rn  (np) 
common /sac  1 /sa ( nsar) 

commori/phiOl/pOl  (np)  /ph i02/p02  ( np )  /ph i03/p03  ( np ) 
c  ommon/ph i 1 1 /pi  1 (np) /phi  12/p  12 (np) /phi 13/p13 (np) 
common/phi21/p21 (np) /phi22/p22 (np) /f hi23/p23 ( np) 
common /phi i l/pi 1 (np)/phii2/pi2<np) 
common/phi i 3/ Pi3 ( np ) /phi i4/pi4 (np ) 
common/phi4/e (np2p) /phi5/f (np2p) 
common /bndv 1 s/bcOp, be Op  1 , be 1p1 , be  Ip, be  1 , bcO 
c ommon/er r c 1 /err  < pp2p ) 
c ommon/ser  rc 1 /amxer 1 , amxer  2, ermax 
common/duml/dl (np2p> /dum2/d2 (np2p) 

commons  for  spline  intesrals 

ssi(np2,i),dsi(np2,7,J),tsi(np2,49,k) 
i=number  of  sinsle  integrals 
J=number  of  double  integrals 
k=number-  of  triple  intesrals 


1976 


c  ommon/ss i c 1 /ss i ( hp2p , nsi p ) 
common/ds i c 1 /dsi (np2p , 7, nd i p ) 
common/tsic  1/tsi  (np2p,  49,  r.tip) 


c  o rrmi o n  / ns s i f / n s i 7  n s j ( n s i f  >  7  ns  v ( 2 7  n s i p ) 
c  ommon  /  n d s  i  f- / n d  i  7  n d J  ( nd  i  p  )  7  n d  v  ( 4 7  n d  i  f  ) 
c  ommon/ nt s i f /nt i 7  nt J ( nt i f ) 7  nt v (67  nt i f  > 


equivalence  (537  sa3) 
dimension  5  a  3  ( 4  7  4  7  n  r  2  f  ) 


: - computes  the  function  values  at  the  test  points  from  the  sF-line  coef 

do  10  i = 1 1 n 

v  <  i)=c  (i)  *f-01  (  i)+c  (  i  +  1 )  *r02(i)+c  <  i+2)  *>03  (  i  > 

10  continue 
return 

entry  repsp(c*v) 
do  11  i  =  1 7  n 

v  <  i)=c  (i  >  *f-11  (i ) +c  (i  +  l>  *F-12(i)+c  (i+2)  *f  13(i> 

11  continue 
return 

entr-Y  r  epspp  (c  7  v) 
do  12  i  =  1 ?  n 

v  (i)=c  (i)  *f-21  <i)+c  <  i  +  1)  #e22<i>  +c  ( i+2)  *f-23  <  i  ) 

12  continue 
retur-n 

entr-Y  rersi  (c  7  v) 


v  < 1 ) =0. 0 
do  13  i=27n 
v(i)=v(i-l)+c(i-l)*Fi 
2  c (i+2) *pi4<i) 

13  continue 
return 
end 

subroutine  rmove(rnew) 
IMPLICIT  REAL #8 ( A-H, 0-Z ) 
*****  written  bY  J.  c.  wileY 


i  3  (  i  )  + 


univ.  of  texas  at  austin  Jan  1976 


PARAMETER (NP 


2 

3 

4 

5 

6 
7 
7 

3 

4 


n  e  9  f-  = 
nsips 
ndir* 
n  t  i  p  3 


217 

1 7 
1 7 
1  7 

4, 


n  2  d  p  =  1 0 1 7 

nm2r=nr-27 nml p=np-l 7  nr  1 f  =np+l 7  of 2f  =  np+27 
ne^s^F  “ne'qF  *ne<qF  7  nep=4*nee»F-+l  7 
maxF-«*ne«»F  *n9dF  7  nsap=16*np2P7  nalF-=7*ne,qsi?F  7 
ndm2F  *riF  2f  7  ndm3ps;ne'qs'qp*nF  2f  7 


ndm23F:=3*ne<HF  7  ndm33F-=3*nessqF  ) 


c  ommon/nc  1  / nm2  7  nml  i  n»  nF- 1 7  hf-2 
c  ommon/ re  1  /  r  <  n  1=  )  ,  r  n  <  hf  ) 
c  0  mmo n/sacl/salnsap) 

c  ommon/F-h  i  01  /f-01  ( nF >  /rh  i  02/f02  <  np )  /  F-h  i03/  f  03  (nr) 
c  orrimon/F-hi  1  1/f1  1  (hf  )  /  f  h  i  12/f  12  (riF- )  /  F  hi  13/f  13  (hf  ) 
c  ommon/F-h  i21  / r-21  (hf  )  /p hi22/p22  ( np  >  /ph  i23/p-23  (  hp  > 
c  ommon/F-h  i  i  1  /  F-i  1  (hf  )  /  F-hi  i2/F-  i2  (rip  ) 
c  ommon  /  p  h  i  1 3/  f-  i  3  ( n f-  )  /  f h  i  i 4/  f-  i  4  ( n f-  ) 
c  ommon/F-h  i4/e  (  i-if-2f  )  /ph  i5/f  ( np2p ) 
c  ommo  n  /  b n  d  v  1  s  /  b  c  Op  7  b  c  Of-  1 7  b  c  1  f-  1 ,  h  c  1  f-  5  b  c  1 7  h  c  0 
c  ommon / er  r  c 1 / er  r ( n  p2p ) 
c  ommon  /  ser-r  c  1  ,'amxer  1 7  smxer-27  er-ma x 
c  ornmon/duml  /d  1  ( np2r  )  /dum2/d2  ( np2p  ) 


c 

c#### 

c#### 

c.#### 

c#### 


commons  for  spline  intesrals 

ss  i  (riF-27  i  )  7  dsi  ( np-27  7>  J  )  7  ts  i  ( nr27  497  k  ) 
i-number  of  sinsle  intesrals 
J “number  of  double  intesrals 
k “number  of  triple  intesrals 


c  ommon/ ss i  c  1  /  ss  i  ( hp2f- 7  nsip) 
c ommon /ds i c 1 /ds i ( np2p 7 7 7  nd i f- ) 
common/tsic  1/tsi  (riF-2F-7497  ntip) 

comrnon/nssip/nsi  7  nsJ  (nsip)  7  nsv  (27  nsip) 
common/ndsip/ndi 7  ndJ (ndip) 7  ndv (47  ndip) 
common/ntsip/nti 7  ntJ (ntip) 7  ntv (67  ntip) 

equivalence  (sa7sa3) 
dimension  sa3 (47 4? np2p) 


dimension  rnew(l) 

do  10  i  =  1 »  n 

r  (  i ) =r  new ( i ) 

r  n  ew ( i ) =0 . 0 . 

continue 

return 

end 

subroutine  sfeval 
IMPLICIT  REAL*8(A-H,0-Z) 


c##*#*  written  bY 


wi  1  eY 


univ.  of  texas  at  austin 


1976 


PARAMETER (NP  = 
2  n  e  s  f  = 


nsiF-=  li 
n  d  i  f  =  1 7 

ntip=  45 
nsdF-  =  101 7 

nm2p=riF--27  nmlp=*rip-l  7  nplp=nr+j ,  np2p=np+2 7 
neqsqp=neqF  fnesp  7  nep=4*neqp  +  l 7 
maxp®neqp#n3dF  7  nsap=16*np2p  7  na lp=7*neqsqp7 
ndm2r“neqp*np2p 7  ndm3p=neqsqF  *np  2f-  7 
ndm23p=3<tneqp 7  ndm33p=3*neqsqp ) 


c  ommon / nc  1  / nm27  n  m  1 7  n  7  n  f  1 7  of  2 
c  ommon /rcl/r  ( n  p )  ,rn(  n  f  ) 
c  o  mm  o  n / s  a  c 1 / s  a  <  n  s  a  p  ) 

c  ommon  / ph  iOl  / f  OI  (nr)  / F-h  i 02/ f  02  ( ri r  )  /  r  h  103/ f-  03  ( np ) 
common/phi  11/f  1 1  (np)  /f  hi  12/f 12  (hf  )  /f  hi  13/f 13  (of  ) 
<•  c  ommon /  p  h  i  2 1  / f-2  1  <  n f  >  /f  h i22/F 22  ( r.F  )  /ph i 23/f 23 < np > 

c  ommon  /  F-h  i  i  1  /  f  i  1  (  riF  )  /  F-h  i  i  2/p  x2  (  of  ) 
common/phi i3/pi3(np) /phi i4/pi4 (np) 
c  ommon/ph i 4/e  ( tip2f  )  /ph i5/f  ( np2p ) 
c  ommon/hnd  v  1  s/bcOp>  be  Op  1  •>  be  1  p  1  ■.  be  1  p  ,  be  1 »  bcO 
c  ommon/erre  1/err  ( riF-2p ) 

0  common/serrc  1/amxerl  ■>  amxer-2>  ermsx 

c  ommon /d urn 1 /dl ( np2p ) /dum2/d2 ( op2p ) 
c 

e  commons  for  spline  integrals 

C 

c####  ssi (np2i i ) t dsi (np2> 7> J ) i tsi (ne2> 49i k ) 

^  c####  i-number  of  single  integrals 

c####  J-number  of  double  integrals 

c####  k=nur.iber  of  triple  integrals 


c.  ommo n/ s s  i  c  1  /  s s  i  <  n p2f  7  ns  i  p ) 
c  ommon/ ds i c 1 /ds i (hp2p  >  7»  ndip) 
c  ommon/ ts i c 1 / 1 s i ( hf2p  7  49»  nt i p ) 

c  ommon/nss i f/ ns i >  nsJ ( ns i f ) »nsv(2»nsip) 
common/ndsip/ndi j  ndj (ndip) >  ndv (4, ndip) 
common/ntsiF/nti , ntJ (ntip) i ntv (6, ntip) 

equivalence  (sa»sa3) 
dimension  sa3 (4, 4, np2p) 


(% 


1* 


c 

c 

dimension  p  1  (5)  »  p2  (5)  »  f-3  (5) 
np3=np2+l 

c - single  spline  integrals55 

do  8  i=l,5 
p2 ( i ) =0. 0 
8  f-3  (  i  )  =0.  0 
f-2  ( 1 )  =1 . 
f-3  ( 1 )  =1 . 

if  (nsi .  eci.O)  go  to  20 
do  10  nsiv=l>nsi 
do  10  i " 1 , n f 2 
lmi n=maxO (5-i > 1 ) 

Imax=min0(rip3-i  »4) 
sum=0.0 

do  11  1  =  1 m i n » 1  ma  x 
idx=4*l+16*i-19 

call  pop ( sa ( i dx) 7  p 1 7  nsv ( 1 7  ns i v) 7  nsv  <2?  ns i v ) 5  ns  1 > 

11  sum=sum+saus  10  (  f-1  *  p2>  p3j  nsJ  (nsi  v)  -ns  1  >  r  (  i  +  1-4)  »  r  (  i  +  1-3) 
10  ss i ( i i ns i v ) -sum 

c - double  spline  integrals. 

20  if(ndi.es.O)  so  to  30 
do  19  i=l,5 
19  f3 ( i ) =0, 0 
f  3  ( 1 )  =  1 . 0 
do  29  nd i v= 1 , nd i 
do  29  i=l,np2 
do  29  vi  =  1  *  7 
i  J  -  i  +J  --4 

d  s  i  ( 3. 7  J  »  n  d  i  v )  =0 . 0 

i f ( i J . 1 t . 1 . or . i j . o t . np2)  go  to  29 


) 


m.\  hi 


mlb»3uJhn*  lE  ;„l 


lmi  n=maxO  (msxO  ( 1  ,  5-i  )  ,  maxO  ( 1  ,  5-i  J  )  +J  -4) 
lmax=minO  (mi nO  (4,  np3-i )  ,  rrnnQ  (4,  np3-iJ  )  +J-4) 
sum=0. 0 

do  21  1=1 min, lmax 

i d  x=4* 1  + 1 6# i - 1 9 

i  d  xk  =4#  (1+4-J  )  +16*  (  i+J-4)  -19 

call  FOF(sa(idx) ,  r  1 ,  n  d  v  ( 1 ,  n  d  i  v ) >  n d  v ( 2 ,  n d i v  >  ,  n  s  1 ) 
call  f  op  (  sa ( i dxk )  ,  F2,ndv(3,ndiv)  , ndv  <4, ndi v)  ,  ns2) 

21  sum=sum+9aus 10  (  pi  ,  r2,  f-3,  ndJ  ( rid  i  v  >  -ns  l-ns2,  r  (  i  + 1-4)  ,  r  (  i  +  1-3)  ) 

29  dsi(i,J,ndiv)=sum 

c - triple  arline  intesrals. 

30  if  (nti.es. 0)  return 
do  31  ntiv=l,nti 

do  31  i-linp2 

do  31  J  =  1 , 7 

do  31  k=l,7 

iJ  =  i+J-4 

i k=i+k-4 

idx=J+(k-l>#7 

tsi  (  i  , idx,ntiv)=0.0 

if  <  iJ  .  It.  1 .  or.  iJ  .  st.  np2)  so  to  31 

i f ( ik . 1 t . 1 . or . ik . st . np2)  so  to  31 

lmin=max0 (maxO ( 1 , 5-i ) , maxO ( 1 , 5-i J ) +J -4, maxO ( 1 , 5-i k ) +k-4) 
lmax=mi nO (mi nO (4 , np3-i > , mi nO ( 4 , np3-i J )  +  J  -4,  mi  nO  <4, np3-i k ) +k -4 ) 
sum=0.0 

if ( lmin. 3t. lmax)  so  to  31 
do  32  1= lmin, lmax 
idxl=4*l+16#i-19 
idx2=4* (1+4-J ) +16*iJ-19 
i  dx3=4* ( 1+4-k ) +16# i k-19 

call  pop(sa(idxl) , p 1 , nt v ( 1 , nt i v) , ntv(2»ntiv>  ,nsl) 
call  ror(sa(idx2) ,p2,ntv(3»ntiv),ntv(4,ntiv)  ,  ns2) 
call  P0p(sa(idx3),p3,ntv(5,ntiv),ntv(6,ntiv),ns3) 

32  sum=sum+sauslO ( p 1 »  p2, p3, nt J (nt i v) -ns l-ns2-ns3, r  < i  +  1-4) , r ( i  +  1-3 
tsi  ( i , i dx , nt i v) =sum 

31  continue 
return 
end 

subroutine  sf- 1  er  r  <  c  •>  er- ) 

C  IMPLICIT  REAL*8(A-H,0-Z) 

c  #  #  #  Kt  %  written  b  y  J.  c.  wiley  univ.  of  texas  at  austin  Jan  1976 
dimension  er(l)»c(l) 

c - relative  error  estimate  of  spline  fit. 

c - note:  the  routine  onlY  returns  a  value  in  er(i)  if 

c -  the  error  is  sreater  than  the  initial  value  of  er-(i). 

c - note:  if  er  is  used  with  routine  srid,  srid  zeroes  er. 

c 

c 

PARAMETER (NP  =21, 

2  nesp=  1 , 

3  nsip=  1, 

4  n  d  i  p  =  1 , 

5  ntir=  4, 

6  ri3dp=101, 

7  nm2p=np-2,  nml p=np-l ,  nF-lp=nr+l ,  np2p=np+2, 

7  nesssp=nesF  *nesp , nep=4#nesp+l , 

3  maxp^nespKinsdp,  nsap=16#np2p ,  na 1 p=7#nessap , 

4  ndm2p=nesp!Hrip2p ,  ndm3p=nesssp#np2p , 

5  ndni23r  =3#nesp,  ndm33r=3#ne«issp ) 
c 

c 

c  ommon/nc  l/nm2,  nml ,  n,  npl ,  riF-2 
c  ommon / r  c  1  / r  ( n f-  )  .  r  n  ( n i-  ) 


c  ommon/  sac  1  -'s a  ( n s a f  ) 

c:  ommc  n/rh  iOl  /pOI  (nr)  /ph i02/pQ2  (  n f  )  /  r  h  i  03/ f  03  (  n f  ) 
c ommon/F h i 1 1/f 1 1 (np) /phi 12/p12 (hf) /phi 13/f 13 (np) 

»  c  ommon/F  h i21  /f-21  ( n f  )  /f  h i 22/ f 22  ( np  )  /ph i23/p23  ( np ) 
c  ommon /rh iil/pil(np)/phii2/pi2(np) 
c ommc  n/ph i  i3/p i3  ( np  )  /Phii4/pi4( np ) 
common  /p  h  i4/e(np2p)/phi5/f(  nF-2p ) 
c ommon/bndv 1  s /be  Of-, be Of  I , be IpI ,  be  If, be  1 ,  bcO 
c ommon/erre  1  /'er  r-  ( n f- 2 p  ) 
c ommon/serr  c  1/amxerl  ,  amxer-2,  ermax 
c  ommon/  duml / d 1  ( np2p )  /  dum2/  d2  (np2p) 
c 

c  commons  for  spline  integrals 

c 

c####  ss  i  ( riF  2,  i )  ,  ds  i  ( riF-2,  7,  J  )  ,  tsi  <  np2>  49,  k  ) 
c####  i =number  of  single  integrals 
c####  J=number  of  double  integrals 
c####  k=number  of  triple  integrals 
c 

c  ommon/ss i c 1 / ss i ( np2p >  nsi p) 
c  ommon/ ds i c 1 /ds i ( np2p ,  7, nd ip) 
common/tsic 1/tsi (np2p, 49,  ntip) 
c 

common/ nss i p/ns i ,  ns J ( ns i f ) , nsv(2,nsip) 
common/ndsip/ndi , ndJ (ndip) , ndv (4, ndip) 
c  ommo  n / n  t  s i p / n  t i , n t J  ( n t  i  p )  ,  n t  v  ( 6 ,  n t  i  p  ) 
c 

equivalence  (sa,sa3) 
dimension  sa3 (4 , 4, nr2p ) 
c 
c 

fcntii 1 , x) =sa3 ( 1 , li i)+x#(sa3(2, 1 , i )  +x# (sa3 (3, 1, i)+x*sa3(4,  1 ,  i) ) ) 
f  va 1=0. 0 
do  11  i  =  1  ■.  n 

fval“fva.l+abs(c(i)#fcn(i,4>r(i))+c(i  +  l)i|cfcn(i  +  l»3»r(i)>  + 

2  c  (  i+2)  *fcn  (i+2, 2,  r  <i)  )  ) 

11  continue 
fval=fval/n 

i  f  ( f  va  1 .  eq .  0.  0)  fval=l. 
do  10  i=2, nml 

er  r  o  r =0 . 02625* ( c <i-l)*<  -sa3 (4, 4, i-1 ) ) 

2  +c(i)  #(sa3(4,4,i)  -sa3(4,3,i>  ) 

3  +c  (i  +  1)  # (sa3(4,3,  i+1 > -sa3 <4, 2, i  +  1) ) 

4  +c < i+2) # (sa3 (4,  2, i+2)-sa3(4, 1, i+2) ) 

5  +c  (i+3)  *(sa3<4,  1,  i+3)  )  > 

6  t (MAX  (r ( i+1 ) -r ( i ) , r < i ) -r ( i-1 ) ) *«3> 
err  or=abs ( error /f va  1 ) 

if (error. at. er ( i ) ) er  < i ) “error 
10  continue 
r  e  t  u  r  n 
end 

subroutine  s p v 1 (i.x.v>m) 

C  IMPLICIT  REAL*8(A-H,0-Z) 

written  bY  J.  c.  wileY  univ.  of  texas  at  austin  Jan  1976 

c - computes  m  values  of  the  i-th  b-spline  at  the  m  x-values  and 

c  returns  them  in  y.  note  that  the  x-values  are  assumed  to  be 

c - ordered. 

c - 

dimension  x(m),Y(m) 
c 
c 

PARAMETER (NP  =  21, 

2  ntip*  1  , 


3  nsip  = 

4  n  d  i  r-  =  It 

5  ntip=  4, 

6  n  s  d  f  -- 1 0 1  t 

7  r,rri2F  =nF  -2t  nrnl  f=tif  -1 ,  hf  1  p=np+l ,  of  2f =hf  +2 t 

7  n  es s q p  =  n es p # n es f-  7  n ep=4# n e q f-  +  1  ? 

3  rri3XF=ne-ctF  #riBdF->  risaF-=16#rip2F  •,  na 1 f =7#neqsqp 7 

4  ridm2p=neqp*rtF  2p  ?  ri  dm3p =n  e q s  q f  * n f  2 f-  ■> 

5  ridm23p=3#neqp ,  ndrn33F  =3*neqsqp ) 
c 

c 

c 

c  omrno n  /  nc  1  / rim2 >  nm  1 »  ri »  n p  1 7  ri f-2 
c o mm  on /r  cl/r  ( n f- )  >  r  n  ( rt p ) 
c.  omnion/ssc  1  /ss  ( n  s  a  f  ) 

c  ommon/ph  iOl  /pOI  ( np  >  /  phi 02/ f-02  ( np )  /ph  i 03/p03  ( np  ) 
c  omrno  n  /  p h  i  1 1  /  p  1 1  ( n p )  /  p h  i  1 2/  p  1 2  ( n p )  /  p h  i  1 3  /  p  1 3  ( n  p  ) 
common/Phi21/p21  (np)  /  f  h  i  22  /  f22  ( r.  r)  /  p  h  i  23  /  f23  ( n  f  > 
c  ommon/ph i i 1 /p i 1 (np) /phi i2/pi2 (np) 
c  omrno  rt  /  p  h  i  i  3/  f  i  3  ( rip )  /  p  h  i  i  4/  p  i  4  ( np ) 
c  ommon/ph .i 4/e  ( np2r- )  /  ph iS/f  ( np2p ) 
c ommort/hridvl  s/ hcOp t  bcOpl  7  be  IpI t  be  1p»  be  1 »  bcO 
common/er re 1/err (np2p) 
common /sen re  1 /amxer 1 , amxer  2, errnax 
c omrno n/duml /d 1 (np2p) Zdum2/d2 < hp2p ) 
c 

c  commons  for  spline  integrals 

c 

c####  ss  i  ( rtF-27  i  )  ,  ds  i  ( nF-27  7 7  J  )  7  ts  i  ( np27  497  k ) 
c####  i=number  of  sin9le  integrals 
c####  J*=riumber  of  double  integrals 
c####  k=n umber  of  triple  integrals 
c 

c omrno n / ss i c 1 /ss i (np2P7  ns i p) 
c  ommon/ds i c 1 /ds i (np2p»7*ndip) 
common/ tsic 1 /ts i (np2P7 49 7  ntip) 
c 

common/nssip/nsi »  nsJ (ns ip) 7  nsv (2t  nsip) 
common/ndsip/ndi 7  ndJ (ndiF  ) 7  ndv (47  ndip) 
cornmon/ntsip/nti  7  ntJ  (ntip)  7  ntv  (67  ntip) 
c 

equivalence  (sai?a3) 
dimension  sa3 (4> 4? np2p) 
c 
c 

f  i3(a3ia27al.a0>:)=ai  ( a  3+z  $  (0«5‘l>a2+zl*  (al/3»0+r^0*25^a0)  )  ) 

f ix (a3, a2, al , aO, z ) =z#z* (0. 5*a3+z * (a2/3. 0+z * (0. 25*a 1+z *0. 2*a0) ) ) 

kk»l 

do  10  J  =  1 7  m 
v  (j  )=0.0 
do  11  k«kk»n 

i  f  ( x  ( J  >  .  1  e .  r  <  k  )  )  go  to  12 

11  continue 

12  l=maxO (2, k ) -i+3 
kk=k 

ifd.lt.l.or.l.gt.4)  so  to  10 
idx=4*l+16*i-16 

v(j)s(sa(  i  dx-3)  +x(j)#(5a(idx-2)+x(J)#(sa(idx-l)+x(J)*sa(idx)  )  )  ) 
10  continue 
return 

entry  s  p v 1 p ( i 7  x  7  v  7  m ) 
k  k  =  1 


Y  ( j  >  =0 . 0 

do  21  k=kh»n 

if (x  <  j ) . le. r  <k) )  so  to  22 

21  continue 

22  l-rr.axO(2,k)-i+3 
kk=k 

if (1. It. 1. or. 1. st. 4)  so  to  20 
idx=4*l+16*i-16 

Y(J)=s.a(idx-2)+x(J)*(2.0*sa(idx-l)+3.0*x(J)*sa(idx)> 
20  continue 
return 

entry  5fvlPF(iiXiY.m) 
k  k  =  1 

do  30  J  =  l)(Ti 

Y  < j  >=0.0 

do  31  k = k k  5  n 

if(x(J).le.r<k>>  so  to  32 

31  continue 

32  l=maxO <2, k  > -i+3 
kk=k 

ifd.lt.l.or.l.st.4)  so  to  30 
idx=4*l+16*i-16 

v  ( J  )  =2.  0*sa  (  i  dx-1 )  +6.  0*sa  ( idx)  *x  ( j  ) 

30  continue 
return 
end 

SUBROUTINE  JOBTIME  (TENMILI. ) 

INTEGER  LISTITEM (3) , TENMILI 
I NTEGER  SYSSGET JP I , STATUS 

LI  ST ITEM  ( 1 )  =  1031*2**16+4 
LISTITEM  (2)  =  '/.LOC  (TENMILI) 

LISTITEM  (3)  =  0  !  7.L0C  (LCPUEL ARSED) 

STATUS3  SYS$GETJP I ( , , , LISTITEM,  ,  ,  ) 

RETURN 

END 

FUNCTION  GET I ME (BUM) 

INTEGER  TENMILI 
CALL  JOBTIME (TENMILI) 

T=FLGAT (TENMILI > / 100. 

GET IME=T 
RETURN 


